A method for identifying the behavior of underwater cracks in a concrete face rockfill dam

Abstract

This invention discloses a method for identifying underwater crack behavior in face-face rockfill dams. The method comprises: obtaining raw monitoring data and as-built drawings of the face-face rockfill dam and performing preprocessing; constructing a permeability-mechanical coupling state model of the rockfill body to obtain a state field and an uncertainty field; analyzing the face-face deformation and stress based on the state field and as-built drawings to obtain the face-face deformation and stress fields; and extracting and visualizing crack features using the deformation and stress fields and underwater detection data to generate a comprehensive crack behavior report. The method introduces a permeability-mechanical coupling state indicator to partition the rockfill body state space, improving the accuracy of deformation field reconstruction and providing a scientific basis for dam safety assessment.

Description

Technical Field

[0001] The invention belongs to the field of water conservancy projects, and in particular to a method for identifying underwater crack properties of a face rockfill dam. Background Art

[0002] Concrete face rockfill dams (CFRDs) are a mainstream form of dam construction today, and their safe operation is crucial to water conservancy projects and the national economy. As a key anti-seepage structure, the concrete face is susceptible to cracking under long-term water pressure and rockfill deformation. Underwater cracks, in particular, are difficult to directly observe, posing a potential threat to dam safety. Accurately identifying the location, morphology, and development trends of underwater cracks is crucial for preventing leakage accidents, extending the dam's service life, and developing scientific maintenance plans.

[0003] Currently, underwater crack detection primarily relies on methods such as sonar scanning, underwater videography, and manual diving. Sonar scanning captures the outline of underwater structures; underwater videography records surface details; and manual diving provides direct observation. In terms of data processing, common methods include image segmentation, contour extraction, and acoustic signal processing to identify crack characteristics. Empirical or semi-empirical models are used to analyze the correlation between deformation and cracks. Furthermore, multi-source data fusion technology is increasingly being applied to underwater structure inspection, improving the comprehensiveness of information acquisition.

[0004] The most critical issue in existing technologies is the lack of a coupled modeling method that can simultaneously consider the interaction between the stress state and permeability of the rockfill. Rockfills exhibit significant differences in permeability under different stress states, particularly in the critical stress region, where nonlinear variations occur. This permeability-mechanical coupling directly affects panel deformation and crack development. Existing studies often separate seepage analysis from stress analysis or use simple linear relationships to characterize their interaction. This fails to accurately characterize the spatiotemporal evolution of the permeability-mechanical relationship of rockfills under complex stress fields, particularly the local nonlinear behavior in regions of drastic stress gradient changes. This technical difficulty severely restricts the accuracy of underwater fracture property identification and urgently requires innovative solutions. Summary of the Invention

[0005] The purpose of the invention is to provide a method for identifying underwater crack properties of a concrete face rockfill dam, in order to solve at least one technical problem existing in the prior art.

[0006] A technical solution, a method for identifying underwater crack properties of a concrete face rockfill dam, comprises the following steps:

[0007] Obtain original monitoring data and as-built drawings of the concrete face rockfill dam, and preprocess them to obtain standardized monitoring data sets and data reliability indicators;

[0008] Based on the standardized monitoring data set, data reliability index and the as-built drawings of the concrete face rockfill dam, the rockfill state field and state uncertainty field are obtained through the rockfill body permeability-mechanical coupling state modeling.

[0009] Using the rockfill state field, state uncertainty field and as-built drawings of the concrete face rockfill dam, the concrete face deformation field and stress field are obtained through concrete face deformation and stress analysis.

[0010] Based on the panel deformation field, panel stress field and pre-stored underwater detection data, a comprehensive crack behavior report is generated through crack feature extraction and visual decision support.

[0011] Beneficial effect: The present invention realizes high-precision modeling of the permeability-mechanical coupling state of the rockfill body, provides a reliable theoretical basis for subsequent face plate deformation and stress analysis and crack feature extraction, and ultimately achieves the purpose of accurately identifying the underwater crack properties of the face plate rockfill dam. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 This is a flowchart of the steps of a method for identifying underwater crack properties of a face rockfill dam provided in an embodiment of the present application.

[0013] Figure 2 This is a flowchart of the steps for modeling the permeability-mechanical coupling state of a rockfill provided in an embodiment of the present application.

[0014] Figure 3 A flowchart of the steps for characterizing the spatial heterogeneity of the rockfill state provided in an embodiment of the present application.

[0015] Figure 4 A flowchart of the steps for obtaining the panel deformation field and panel stress field through panel deformation and stress analysis provided in an embodiment of the present application. DETAILED DESCRIPTION

[0016] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0017] It should be noted that to clearly illustrate the steps of this application, serial numbers are assigned to each step in the specification. These serial numbers are for illustrative purposes only and do not limit the order in which the steps must be executed. In actual operation, depending on the technical requirements of the specific implementation scenario, the steps may be executed in a different order than shown in the specification, and in some cases, parallel processing between steps may be implemented.

[0018] like Figure 1 As shown, a method for identifying underwater crack properties of a concrete face rockfill dam includes the following steps:

[0019] S1. Obtain original monitoring data and as-built drawings of the concrete face rockfill dam, and pre-process them to obtain standardized monitoring data sets and data reliability indicators;

[0020] Specifically, the original monitoring data includes: structural geometry data: 3D geometric information of the above-water section of the slab, the dam crest, and the downstream dam surface, acquired through 3D laser scanning; underwater sonar scanning data: information on the underwater dimensions and blanket configuration of the slab; deformation monitoring data: information on dam deformation acquired through water-tube sedimentation meters; mechanical parameter data: including stress and strain measurements; hydraulic parameter data: including hydraulic monitoring data such as pore water pressure and flow velocity; temperature field data: information on the temperature distribution inside and outside the dam; and time series monitoring data: a historical record of the changes in various parameters over time. The as-built drawings of the CFRD are the final design and construction records after construction is completed. These drawings detail the design parameters of the dam structure and its components.

[0021] S2. Based on the standardized monitoring data set, data reliability index and the as-built drawings of the concrete face rockfill dam, the rockfill state field and state uncertainty field are obtained through the rockfill body permeability-mechanical coupling state modeling;

[0022] Specifically, the model simultaneously considers water infiltration and the mechanical behavior of the dam structure (rockfill). In other words, the flow of water within the dam (infiltration) and the deformation of the dam under pressure (mechanical behavior) are coupled and influence each other. Numerical simulations yield an "actual dam state field," encompassing information such as stress, deformation, and internal water pressure in each region, creating a "map" of the dam's internal conditions. Because both monitoring data and model parameters have a degree of uncertainty, the model also assesses the reliability of the prediction results for each region. This creates a state uncertainty field, which reveals where the assessment results are highly reliable and where more caution may be required due to data or parameter issues.

[0023] S3. Using the rockfill state field, state uncertainty field, and as-built drawings of the concrete face rockfill dam, obtain the concrete face deformation field and concrete face stress field through concrete face deformation and stress analysis;

[0024] Specifically, numerical simulation or theoretical calculation is performed using the state information of the rockfill as input conditions, combined with the geometric shape and material properties of the components, to calculate the deformation (such as bending, displacement) and stress distribution (such as compressive stress, tensile stress) of the dam panel under the influence of internal driving force, reaction and external environment.

[0025] S4. Based on the panel deformation field, panel stress field and pre-stored underwater detection data, a comprehensive report on crack properties is generated through crack feature extraction and visual decision support.

[0026] Specifically, crack characteristics are extracted and their development mechanisms analyzed, generating a multi-dimensional visualization of crack behavior and providing targeted repair recommendations. In other words, the panel deformation and stress fields can be used to identify areas of abnormality. For example, a sudden increase in deformation or stress concentration in a particular area may indicate the formation or expansion of cracks in a certain part of the structure. Comparing or combining these abnormal areas with underwater inspection data can further confirm and refine the crack characteristics (such as location, length, width, and depth). After extracting crack-related features, data visualization techniques are used to intuitively present the distribution and severity of the cracks through charts, maps, or 3D models.

[0027] This embodiment achieves comprehensive and dynamic monitoring of the dam's structural safety status through multi-source data integration, precise physical coupling modeling, and intuitive visual decision support. From reducing data errors and improving model accuracy to timely and scientific identification of potential structural problems, the entire process has enhanced the technical level of dam health monitoring and risk warning.

[0028] According to one aspect of the present application, the pre-processing steps are specifically as follows:

[0029] S11, reading 3D laser scanning data, underwater sonar scanning data, and water tube sedimentation meter monitoring data, and calibrating them through time tags to generate original monitoring data;

[0030] S12. Read the original monitoring data, identify and mark outliers using statistical methods, and calculate the data quality index of each monitoring point;

[0031] S13, read the original monitoring data and data quality indicators, perform spatiotemporal coordinate conversion and alignment of multi-source data, and generate a spatiotemporal aligned dataset;

[0032] S14. Read the spatiotemporal aligned dataset, apply the adaptive wavelet noise reduction algorithm according to the characteristics of different sensors and environmental interference factors, and output the standardized monitoring dataset and the data reliability index for comprehensive evaluation.

[0033] like Figure 2 As shown, according to one aspect of the present application, the steps of obtaining the rockfill body state field and state uncertainty field by rockfill body permeability-mechanical coupling state modeling include:

[0034] S21. Based on the as-built drawings of the concrete face rockfill dam and the structural information in the standardized monitoring data set, a physical model of the rockfill body is established, and a geometric model of the rockfill body and material partition information are generated;

[0035] S22. Combining standardized monitoring data sets, data reliability indicators, and material partitioning information, characterize the spatial heterogeneity of the rockfill state and form a spatial distribution feature of the rockfill state;

[0036] S23. Based on the spatial distribution characteristics of the rockfill state and the time series data in the standardized monitoring data set, simulate the dynamic evolution of the time-varying state and construct a spatiotemporal evolution model of the rockfill state;

[0037] S24. Integrate the spatiotemporal evolution model of the rockfill state and the data reliability index, analyze and quantify the uncertainty and its propagation law, and output the rockfill state field and state uncertainty field.

[0038] Specifically, the system retrieves standardized monitoring datasets, data reliability indicators, and material partitioning information, and employs a stress-state-sensitive self-organizing mapping method to identify and characterize the spatial heterogeneity of the rockfill's permeability-mechanical coupling characteristics, outputting the spatial distribution characteristics of the rockfill state. Based on the spatial distribution characteristics of the rockfill state and the time series data from the standardized monitoring dataset, a time-varying evolution model of the rockfill's permeability-mechanical coupling state is constructed, outputting a spatiotemporal evolution model of the rockfill state. Based on the spatiotemporal evolution model of the rockfill state and the data reliability indicators, a Monte Carlo simulation and variational Bayesian inference method are used to quantify the uncertainty of the rockfill state parameters and their propagation patterns, generating a rockfill state field and a state uncertainty field.

[0039] This embodiment extracts fine structural information from static design drawings and real-time monitoring data to construct a geometric model and material partitioning that truly reflects the complex internal structure of the dam, ensuring that subsequent calculations are based on a solid and accurate information foundation; by characterizing the heterogeneity of the rockfill state space, it reveals local differences in internal materials, structures, and stress conditions, enhancing the ability to locate risks and manage regional differentiation; the model can capture the dynamic evolution trend of the dam state over time, providing early warning and trend analysis support for fault prediction, accident prevention, and maintenance decision-making; it not only provides a prediction of the current state, but also clarifies the uncertainty range of the prediction results, scientifically supporting engineers in focusing on monitoring risk areas and formulating response measures; it realizes all-round monitoring and description of the dam rockfill from spatial distribution to temporal evolution.

[0040] like Figure 3 As shown, according to one aspect of the present application, the steps of characterizing the spatial heterogeneity of the rockfill state and forming the spatial distribution characteristics of the rockfill state include:

[0041] Extract relevant monitoring data from the standardized monitoring dataset and construct a permeability-mechanical coupling state index that considers effective stress, pore water pressure, and strain state;

[0042] Based on the permeability-mechanical coupling state index and material partition information, stress state sensitive self-organizing mapping is performed to divide the rockfill state characteristic zones;

[0043] Based on the zoning of rockfill state characteristics, a regional mechanical characteristic dynamic model considering stress-seepage coupling is established to generate a set of local mechanical response models;

[0044] Based on the local mechanical response model set and data reliability indicators, boundary adaptive fusion is implemented to achieve seamless connection of different partition models and output the spatial distribution characteristics of the rockfill state.

[0045] Specifically, the local mechanical response model set and data reliability index are read, and the boundary adaptive fusion method is used to achieve seamless connection of different partition models, and output the spatial distribution characteristics of the rockfill state. The boundary adaptive fusion algorithm is: CSI(x) = Σi[Wi(x)·CSIi(x)]; Wi(x) = Ri·ρi(x) / Σj[Rj·ρj(x)]; ρi(x) = exp(-di(x) 2 / 2li 2 ); CSI(x) is the fused state indicator field; CSIi(x) is the state indicator of the i-th local model; Wi(x) is the weight of the i-th model at position x; Ri is the reliability score of the i-th model; ρi(x) is the distance decay function; di(x) is the distance from point x to the center of the i-th region; li is the characteristic length scale of the i-th region; x is the spatial coordinate point; i and j are the indices of the state characteristic partitions. A boundary-adaptive fusion method is used to achieve seamless connection between models in different regions. Local models are weighted and fused based on data reliability and spatial distance, effectively addressing the spatial heterogeneity of the permeability-mechanical relationship within the rockfill.

[0046] This embodiment comprehensively analyzes the mutual influence of stress, seepage, and deformation states of the rockfill body, accurately characterizing its spatial heterogeneity; dynamically simulates the response characteristics of the rockfill body under different stress states, thereby improving prediction and control capabilities; connects different partition models to ensure seamless connection between partitions, which helps to achieve overall optimization and coordination of the spatial distribution characteristics of the rockfill body state; and enhances the credibility and accuracy of the technical results during the model generation process, ensuring that the final state distribution characteristics can accurately reflect the actual situation.

[0047] Because traditional methods typically only consider a single indicator, water content or pore water pressure, to characterize the degree of wetting, ignoring the complex coexistence of solid, liquid, and gas phases within the rockfill, according to one aspect of the present application, the steps for constructing a permeability-mechanical coupling state indicator include:

[0048] Extract effective stress, pore water pressure, and strain state data from standardized monitoring datasets;

[0049] Establish a mathematical model that comprehensively considers effective stress, pore water pressure and strain state, and calculate the permeability-mechanical coupling state index CSI(x, t);

[0050] CSI(x,t) = α·σ'(x,t) / σ'max + β·P(x,t) / Pmax + γ·ε(x,t) / εmax;

[0051] Where CSI(x, t) represents the permeability-mechanical coupling state index of the spatial point x at time t; σ'(x, t) represents the effective stress; P(x, t) represents the pore water pressure; ε(x, t) represents the strain; σ'max represents the maximum effective stress; Pmax represents the maximum pore water pressure; εmax represents the maximum strain; α, β, and γ represent weight coefficients, and α+β+γ=1; x represents the spatial coordinate point; and t represents the time point.

[0052] This embodiment introduces a comprehensive seepage-mechanical coupling state index (CSI). Using the formula CSI(x, t) = α·σ'(x, t) / σ'max + β·P(x, t) / Pmax + γ·ε(x, t) / εmax, the three key parameters of effective stress, pore water pressure, and strain are unified into a quantitative framework. This allows for a comprehensive characterization of the complex internal state of the rockfill, overcomes the limitation of traditional methods that separate seepage and mechanical analysis, and addresses the problem that a single indicator cannot fully describe the wetting state.

[0053] The traditional SOM algorithm cannot effectively identify the critical point of phase transition of the wetting degree when processing the wetting zoning of the rockfill body. Therefore, according to one aspect of the present application, the steps of performing stress state sensitive self-organizing mapping to divide the rockfill body into state characteristic zones include:

[0054] Input the permeability-mechanical coupling state index and material partition information;

[0055] Applying the improved self-organizing map algorithm to partition the rockfill state space and generate rockfill state characteristic partitions;

[0056] The improved self-organizing map algorithm contains the following core update formula:

[0057] wij(t+1) = wij(t) + η(t)·h(i,BMU(x),t)·[x - wij(t)]·ψ(σ'(x,t));

[0058] where wij(t) is the j-th dimension weight of the i-th neuron at time t; η(t) is the learning rate, η(t) = η0·exp(-t / τ1), η0 is the initial learning rate; h(i, BMU(x), t) is the neighborhood function, h(i, BMU(x), t) = exp(-||ri -rBMU(x)|| 2 / (2σ 2 (t))); ri is the position of the i-th neuron in the grid; r is the spatial position; σ(t) is the neighborhood radius, σ(t) = σ0 exp(-t / τ2); τ1 and τ2 are time constants; σ0 is the initial neighborhood radius; BMU(x) is the best matching unit index; ψ(σ'(x, t)) is the stress state sensitivity function, ψ(σ') = [1 + exp(-k (σ'-σ'0))] -1 ; σ'(x, t) is the effective stress; σ'0 is the stress critical value; k is the stress sensitivity parameter; x is the eigenvector of the input sample, rBMU(x) is the position of the best matching unit in the grid, and t is the number of iterations.

[0059] The stress state sensitivity function ψ(σ') is used to adjust the sensitivity of the algorithm to different stress state areas and improve the recognition accuracy of critical stress areas.

[0060] This embodiment proposes a stress state sensitive self-organizing map algorithm, which introduces a stress sensitive function ψ(σ') = [1 + exp(-k·(σ'-σ'0))] into the traditional SOM algorithm. -1 , achieving accurate identification of critical stress regions. This enables the algorithm to have differentiated response sensitivities to regions with different stress states, particularly enhancing the ability to identify regions with dramatic stress gradient changes, solving the problem that traditional methods have difficulty accurately capturing nonlinear change characteristics.

[0061] Traditional wetting diffusion models typically use homogeneous assumptions and fixed diffusion coefficients, which cannot capture the nonlinear characteristics of the wetting process within the rockfill. Therefore, according to one aspect of this application, the steps of establishing a regional mechanical property dynamic model include:

[0062] Analyze the mechanical characteristics of the rockfill state characteristic zones;

[0063] A dynamic model of regional mechanical characteristics considering stress-seepage coupling is constructed: ∏CSI(x, t) / ∏t = ▽·[D(x, σ', CSI)·▽CSI(x, t)] + Q(x, t, σ');

[0064] where ∏ is the partial derivative; D(x, σ', CSI) represents the diffusion coefficient tensor, D(x, σ', CSI) = D0 exp(-α eσ')·exp(λ·CSI)·R(θ); R(θ) represents the anisotropic rotation matrix; D0 represents the reference diffusion coefficient; α e represents the effective stress influence coefficient; λ represents the coupling state sensitivity coefficient; Q(x, t, σ') represents the source-sink term, which represents the state change rate caused by local stress change; θ represents the angle between the main diffusion direction and the reference direction; CSI(x, t) is the permeability-mechanical coupling state index of the spatial point x at time t; σ' is the effective stress; ▽ is the gradient operator; x is the spatial coordinate point; t is the time point; ▽· is the divergence operator.

[0065] This embodiment constructs a regional mechanical property dynamics model, in which the diffusion coefficient tensor simultaneously considers the influence of effective stress and coupling state indicators. It can accurately describe the nonlinear mechanical response of the rockfill under complex stress states and solves the problem of insufficient accuracy caused by the use of a fixed diffusion coefficient in traditional diffusion models.

[0066] According to one aspect of this application,

[0067] The steps to generate a comprehensive report on fracture properties through fracture feature extraction and visual decision support include:

[0068] Read the panel deformation field and panel stress field, calculate the crack sensitivity index of each area of the panel based on crack mechanics theory, and obtain a crack sensitive area map;

[0069] Read underwater detection data and crack sensitive area maps, apply wavelet transform domain morphological decomposition and adaptive reconstruction methods to identify crack location, direction and width, and obtain crack characteristic data;

[0070] Read crack characteristic data, panel stress field, and panel deformation field, comprehensively analyze the cause, development status, and risk level of the cracks, and obtain the crack property determination results;

[0071] Read the crack behavior determination results, crack characteristic data and crack sensitive area map, generate a multi-dimensional crack behavior visualization expression, provide targeted maintenance suggestions based on risk analysis, and output a comprehensive crack behavior report.

[0072] Among them, the wavelet transform domain morphological decomposition and adaptive reconstruction method includes:

[0073] Perform local illumination normalization and color channel weighted fusion on the underwater detection data to obtain a fused grayscale image;

[0074] Apply multi-layer wavelet decomposition to the fused grayscale image to obtain low-frequency approximation coefficients and multi-layer detail coefficients;

[0075] Calculate the directional coherence measure and phase consistency of the detail coefficients to obtain a directional coherence measure set and a phase consistency map;

[0076] An adaptive enhancement function is constructed based on the directional coherence measure set and the phase consistency map to perform selective enhancement on the wavelet coefficients to obtain enhanced wavelet coefficients.

[0077] Perform inverse wavelet transform and local statistical enhancement on the enhanced wavelet coefficients to obtain an enhanced crack image;

[0078] Morphological sharpening and crack feature measurement are applied to the enhanced crack image to obtain crack feature data.

[0079] The calculation of the directional coherence measure includes: reading the horizontal detail coefficient and the vertical detail coefficient, and calculating the initial directional coherence measure for each scale j and spatial position (x,y): M_j(x,y) = |∑ω cH_j(x,y)·cH_j(x+ω,y) + cV_j(x,y)·cV_j(x,y+ω)| / sqrt(cH_j(x,y) 2 + cV_j(x,y) 2 ); calculate the directional consistency compensation coefficient in the local area: γ_j(x,y) = max(0, 1 - var(θ(x+Δx,y+Δy)) / π), where θ represents the local gradient direction and var represents the variance; combine the initial directional coherence measure and the directional consistency compensation coefficient to calculate the compensated directional coherence measure: M'_j(x,y) = M_j(x,y)·(1 + α_j·γ_j(x,y)), where α_j is the scale adaptation coefficient; ω is the neighborhood range; cH_j(x,y) represents the horizontal detail coefficient of the spatial position (x, y) at scale j; cV_j(x,y) represents the vertical detail coefficient of the spatial position (x, y) at scale j; Δx,y is the local offset of the spatial position (x, y).

[0080] According to one aspect of the present application, the construction of the adaptive enhancement function includes:

[0081] Read the directional coherence measure set and phase consistency map, and construct an adaptive enhancement function: f_j(M',PC') =(1 + tanh(α_j·(M'_j - β_j))) / 2·(1 + γ_j·PC'); adaptively determine the enhancement function parameters: α_j =α_0·(1 + Δ·SNR_j); β_j = β_0·(1 - ρ·SNR_j); γ_j = γ_0·(1 + ω·PC_mean), where SNR_j is the estimated j-th scale signal-to-noise ratio, PC_mean is the average phase consistency, β_j is the threshold parameter at the j-th scale; PC' is the local phase consistency value; α_0 is the initial scale adaptation coefficient; β_0 is the baseline value of the initial directional coherence measure; ρ is the adjustment coefficient related to the signal-to-noise ratio; γ_0 is the initial phase consistency compensation coefficient. Apply the enhancement function to the wavelet detail coefficients: cH'_j(x,y) = cH_j(x,y)·f_j(M'_j(x,y), PC'(x,y)); apply the amplitude limit function to prevent over-enhancement: cH''_j = sign(cH'_j)·min(|cH'_j|, λ_j·std(cH_j)), where λ_j is the scale adaptation coefficient.

[0082] According to one aspect of the present application, step S23 may also be:

[0083] S231. Multi-scale state decomposition and transition point identification.

[0084] The spatial distribution characteristics of the rockfill state, S_space, and time series data, T_data, from the standardized monitoring dataset are read. Wavelet decomposition is then applied to the state time series, CSI(x,t), at each spatial point x. The Daubechies wavelet (db4) is used as the basis function, and a five-layer discrete wavelet transform is performed: CSI(x,t) = ∑j∑k c_j,k(x)ψ_j,k(t) + ∑k a_J,k(x)φ_J,k(t); where j = 1, 2, ..., J represents the scale level (J = 5), k is the time displacement, ψ_j,k(t) is the wavelet function, φ_J,k(t) is the scaling function, c_j,k(x) is the detail coefficient, and a_J,k(x) is the approximation coefficient. The output is a multiscale wavelet coefficient set, W = {a_J,k(x), c_j,k(x)}.

[0085] Read the multi-scale wavelet coefficient set W and calculate the energy distribution of each scale wavelet coefficient E_j(x,t) = ∑k |c_j,k(x)| 2 (within a window around time t); E_J(x,t) = ∑k |a_J,k(x)| 2 (within the window around time t);

[0086] For each scale j, the energy mutation rate R_j(x,t) = |E_j(x,t) - E_j(x,t-Δt)| / max(E_j(x,t-Δt), ε) is calculated, where Δt is the time step and ε is a small constant to prevent division by zero. To enhance meaningful mutations, an energy density modulation R'_j(x,t) = R_j(x,t)·(1 + Δ·log(E_j(x,t) / E_mean)) is introduced, where E_mean is the mean energy and Δ=0.3 is the modulation coefficient. The modulated energy mutation rate set {R'_j(x,t)} is output. The modulated energy mutation rate set {R'_j(x,t)} and the physical parameters P_data from the standardized monitoring dataset are read to extract key physical quantities: effective stress σ'(x,t) and pore water pressure p(x,t). Calculate the rate of change of physical parameters: Rσ(x,t) = |σ'(x,t) - σ'(x,t-Δt)| / σ'_max; Rp(x,t) = |p(x,t) - p(x,t-Δt)| / p_max; construct the physical trigger correlation function C_j(x,t) = max(corr(R'_j(x,t±τ), Rσ(x,t)), corr(R'_j(x,t±τ), Rp(x,t))). Where corr represents the sliding window correlation coefficient, and τ is the time lag / lead range (±10 time steps). This function captures the time lag between state changes and physical trigger conditions. Output the physical trigger correlation set {C_j(x,t)}.

[0087] The set of modulation energy mutation rates {R'_j(x,t)} and the set of physical trigger correlations {C_j(x,t)} are read, and an adaptive threshold function τ_j(x,t) = μ_j(x) + α·σ_j(x,t)·(1 + β·exp(-γ·t_s(x)))·(1 -ω·C_j(x,t)) is designed. Here, μ_j(x) is the base threshold (local average mutation rate), σ_j(x,t) is the local standard deviation, t_s(x) is the time since the last detected transition point, and α=2.5, β=1.0, γ=0.1, and ω=0.3 are control parameters. A physical correlation modulation term (1-ω·C_j(x,t)) is introduced to make mutations with high physical correlation more likely to be detected as transition points. Identify the set of transition points TP_j(x) = {t | R'_j(x,t) > τ_j(x,t) and local maxima}. Perform multi-scale transition point fusion to eliminate redundancy and false positives: TP(x) = Filter({TP_j(x)}, d_min, n_min). Filter filters transition points based on the minimum time interval d_min and the minimum number of supported scales n_min to retain the most significant true transition points. Output the state transition point set TP(x) and the transition point feature set T_char.

[0088] S232. Multi-scale hierarchical physical modeling.

[0089] The multiscale wavelet coefficient set W and the rockfill state spatial distribution characteristics S_space are read to extract the low-frequency approximate coefficients a_J,k(x). A physical diffusion model is constructed: da_J(x,t) / dt = ▽·[D_L(x,φ(x,t))·▽a_J(x,t)] + S_L(φ(x,t)), where d is the sign of the partial derivative; D_L is the low-frequency diffusion coefficient tensor, φ(x,t) is the physical state vector (including stress, pore pressure, etc.), and S_L is the source and sink term. The diffusion coefficient D_L(x,φ) = D_0·exp(-α_e·σ'(x,t))·exp(λ·a_J(x,t))·R(θ(x)); where D_0 is the reference diffusion coefficient, α_e is the effective stress influence coefficient, λ is the coupled state sensitivity coefficient, R(θ) is the anisotropic rotation matrix, and θ is the principal diffusion direction. The source-sink term S_L represents the equilibrium state drive, S_L(φ) = r_L·(a_J,eq(φ) - a_J); where r_L is the relaxation rate and a_J,eq(φ) is the equilibrium value for a given physical state φ. A finite difference method is used for discrete solution, outputting the low-frequency state evolution model M_L.

[0090] Read the multi-scale wavelet coefficient set W, the state transition point set TP(x), and the transition point characteristic set T_char, and construct a scale adaptive model dc_j(x,t) / dt = D_j·▽ for each detail scale j 2 c_j(x,t) + v_j ▽ c_j(x,t) + S_j(φ,a_J,t); the source term S_j is decomposed into three parts: S_j(φ,a_J,t) = S_j,eq(φ,a_J) + S_j,trans(t) + S_j,noise; the equilibrium term S_j,eq represents the drive towards the equilibrium state, S_j,eq(φ,a_J) = r_j · (c_j,eq(φ,a_J) - c_j); the transformation term S_j,trans is activated near the detected transition point, S_j,trans(t) = ∑_{t_c∈TP} A_j(t_c) · h_j(t-t_c); where h_j is the transformation function, h_j(t) = (1-exp(-ρ_j·t)) / (1+exp(-ρ_j·t)); A_j(t_c) is the transition amplitude, which is related to the strength of the physical trigger condition. The noise term S_j,noise is a random perturbation that simulates small-scale random fluctuations. Output is the multi-scale detail evolution model set {M_j}.

[0091] The low-frequency state evolution model M_L and the multiscale detail evolution model set {M_j} are read to construct a cross-scale coupling mechanism. First, the inter-scale influence matrix I_{i,j} = influence intensity (from scale i to scale j) is established; the influence intensity is estimated using wavelet partial causality analysis based on historical data. Second, the cross-scale feedback term F_j(t) = ∑_i I_{i,j}·∫ K_{i,j}(t-τ)·c_i(τ)dτ is designed, where K_{i,j} is the inter-scale influence kernel function, representing the time delay effect. Finally, the evolution model dc_j(x,t) / dt = ... + F_j(t) is updated. Simultaneously, the equilibrium value a_J,eq(φ) in the low-frequency model is modified to a_J,eq(φ,{c_j}) to make it dependent on the detail coefficients, forming a bidirectional coupling. The cross-scale coupling coefficient set C_cross and the updated multi-scale evolution model set M = {M_L, {M_j}} are output.

[0092] S233. Physical constraint reconstruction and uncertainty quantification.

[0093] Read the multi-scale evolution model set M and the deformation physical constraint condition set C_phys, and construct the physical consistency constraint L[CSI] - f(φ) = 0; where L is the mechanical differential operator and f(φ) is the body force term. Transform it into the optimization problem min_{a',c'} ||a'-a|| 2+ ∑_j||c'_j-c_j|| 2 + λ||L[CSI'] - f(φ)|| 2 where a' and c'_j are the optimized wavelet coefficients, and CSI' is the corresponding reconstructed state field. The constraint weight λ(x) = λ_0·(1 + κ·Ξ(x)); where Ξ(x) is a physical importance index, taking a high value in areas of stress concentration and large pore pressure gradients. The optimization problem is solved using the alternating direction method of multipliers (ADMM), outputting the physically consistent optimization coefficients {a'_opt, c'_j, opt}.

[0094] Read the physical consistency optimization coefficients {a'_opt, c'_j,opt} and the multi-scale evolution model set M to perform state evolution prediction. Introduce Bayesian posterior analysis p(θ|D) ∝ p(D|θ)p(θ); where θ is the model parameter vector and D is the historical observation data. Use the Markov Chain Monte Carlo (MCMC) method to sample the posterior distribution and generate the parameter sample set {θ i}. Execute model prediction for each sample to obtain the state evolution sample set {CSI i (x,t)}. Calculate the predicted mean and variance: CSI_mean(x,t) = (1 / N)·∑_i CSI i (x,t);CSI_var(x,t) = (1 / N)·∑_i (CSI i (x,t) - CSI_mean(x,t)) 2 ; Output state prediction field CSI_pred and prediction uncertainty field U_pred.

[0095] Read the state prediction field CSI_pred, the prediction uncertainty field U_pred, and the time series data T_data from the standardized monitoring dataset to perform sensitivity analysis and model validation. First, calculate the relative sensitivity index S_θj = (θj / CSI)·(dCSI / dθj) to determine the degree of influence of key parameters on the prediction results. Second, perform historical data validation RMSE = sqrt((1 / N)·∑ (CSI_pred - CSI_obs) 2 );R 2 = 1 - ∑(CSI_pred - CSI_obs) 2 / ∑(CSI_obs- mean(CSI_obs)) 2; Finally, calculate the prediction reliability index RI(x,t) = (1 - U_pred(x,t) / CSI_pred(x,t))·(1 - RMSE / CSI_range); output the model validation report R_validation and the prediction reliability field RI.

[0096] S234, spatiotemporal evolution model integration.

[0097] Read the physical consistency optimization coefficients {a'_opt, c'_j,opt} and the multiscale evolution model set M, and perform multiscale spatiotemporal reconstruction. First, initialize the current state: a_cur = a'_opt; c_j,cur = c'_j,opt; then, step by step according to the evolution model, with time step Δt: a_new = a_cur + Δt·F_L(a_cur, {c_j,cur}, φ); c_j,new = c_j,cur + Δt·F_j(a_cur, {c_i,cur}, φ); where F_L and F_j are the corresponding evolution functions. At each time step, physical constraint adjustment is performed: {a_new, {c_j,new}} = PhysicalAdjust(a_new,{c_j,new}, C_phys); finally, the state field is reconstructed through inverse wavelet transform: CSI(x,t) = IDWT(a_new, {c_j,new}); and the spatiotemporal evolution prediction field CSI_evol is output.

[0098] Read the spatiotemporal evolution prediction field CSI_evol and the set of state transition points TP(x) and perform dynamic reshaping of state partitions. The impact of state transitions on partitions is introduced. First, the local transition saliency index TS(x) = ∑_{t_c∈TP(x)} w(t-t_c)·Magnitude(t_c) is calculated, where w is the temporal weight function and Magnitude is the transition strength. Next, the state partition boundaries are adjusted using an adaptive level set method: dΓ / dt = V_n·N, where Γ represents the partition boundary and V_n is the normal velocity, which depends on TS(x) and the state gradient. Dynamic state partitions D_zones are output.

[0099] Read the spatiotemporal evolution prediction field CSI_evol, the prediction uncertainty field U_pred, and the dynamic state partitions D_zones to generate the final rockfill state field and uncertainty field. Perform partition boundary smoothing CSI_smooth(x,t) = ∑_zw_z(x)·CSI_z(x,t); where w_z(x) is the weight of partition z at position x, to achieve smooth transition in the partition boundary area. Calculate the comprehensive uncertainty field U_total(x,t) = sqrt(U_pred(x,t) 2 + U_model(x,t) 2 + U_bound(x,t) 2 ); where U_model is the model structure uncertainty and U_bound is the partition boundary uncertainty. Output the rockfill state field CSI_field and the state uncertainty field U_field.

[0100] A framework for separate modeling of low-frequency components and detail coefficients was established, and a cross-scale coupling mechanism was used to organically integrate dynamic characteristics at different time scales. This approach is particularly well-suited to the simultaneous slow trends and rapid transitions in the state evolution of rockfill bodies. By constructing a physically consistent optimization problem for wavelet coefficients, a deep fusion of the state evolution model and physical laws was achieved, ensuring that predictions not only conform to historical data trends but also meet physical constraints.

[0101] like Figure 4 As shown, according to one aspect of the present application, the steps of obtaining the panel deformation field and the panel stress field through panel deformation and stress analysis include:

[0102] S31, reading deformation monitoring data and rockfill state field in the standardized monitoring data set, reconstructing discrete monitoring point data using a multi-scale physical constraint kernel function interpolation method, and generating a continuous distribution of the deformation field;

[0103] S32. Based on the rockfill state field, continuous distribution of deformation field, and as-built drawings of the CFRD, a CFRD-rockfill contact mechanics model considering the influence of wetting is constructed to calculate the contact surface stress distribution and output the contact stress distribution.

[0104] S33. Based on the contact stress distribution, the continuous distribution of the deformation field, and the panel structure information in the as-built drawings of the face rockfill dam, analyze the stress state of the panel under complex loads and calculate the panel stress field;

[0105] S34. Integrate the continuous distribution of the rockfill state field, state uncertainty field and historical deformation field, establish a deformation prediction model and verify it with the measured data, and output a verified panel deformation field and prediction accuracy evaluation report.

[0106] This embodiment can generate a high-precision continuous distribution of the deformation field, improving data integrity and spatial resolution; more realistically reflect the dynamic changes in the stress distribution of the contact surface, thereby enhancing the ability to analyze the structural behavior of the panel under the action of wetting; comprehensively analyze the stress state under complex loads, accurately obtain the panel stress field, optimize structural design and construction planning; improve prediction accuracy, and provide reliable judgment on the dynamic change trend of panel deformation.

[0107] According to one aspect of the present application, the steps of reconstructing discrete monitoring point data by a multi-scale physical constraint kernel function interpolation method to generate a continuous distribution of the deformation field include:

[0108] Analyze the as-built drawings of the concrete face rockfill dam and the state field of the rockfill body, extract the deformation physical constraints based on the principle of mechanical equilibrium, and form a set of deformation physical constraints;

[0109] Perform variational mode decomposition on the deformation monitoring data in the standardized monitoring dataset to identify the multi-level characteristic scales of the deformation field and generate a multi-scale characteristic spectrum;

[0110] Combining the deformation physical constraint condition set and the multi-scale characteristic spectrum, a multi-scale kernel function that satisfies the physical constraint is constructed, and a physical response kernel function set is established;

[0111] The physical response kernel function set is used to process the deformation monitoring data in the standardized monitoring data set, and multi-scale kernel function interpolation calculation is performed to achieve high-precision reconstruction of discrete monitoring point data into a continuous deformation field, and output a high-precision continuous distribution of the deformation field.

[0112] Among them, the physical response kernel function set ensures that the interpolation results meet both data fitting requirements and physical constraints, effectively solving the shortcomings of traditional pure mathematical interpolation methods in terms of physical rationality.

[0113] Specifically, the as-built drawings of the concrete face rockfill dam and the state field of the rockfill body are read, and the physical constraints of the deformation field are extracted based on the principle of mechanical equilibrium to generate a set of physical constraints for deformation. The mathematical expression of the physical constraints is: L[u(x)] = f(x, MI(x)); B[u(x)] = g(x) on ∏Ω; L is the differential operator of elastic mechanics, L = μ▽ 2 (·) + (λ+μ)▽(▽·(·)); u(x) is the displacement field; f(x, MI(x)) is the body force, which is related to the wetting index MI(x); B is the boundary operator; g(x) is the boundary condition; λ and μ are the Lame constants, which are related to the material properties; ∏Ω is the boundary of the computational domain; ▽ is the gradient operator; ▽· is the divergence operator; ▽ 2 is the Laplace operator; x is the spatial coordinate point, and on represents the scope of application of the mathematical conditions or the area of action of the constraints.

[0114] This embodiment successfully converts discrete monitoring point data into a high-precision continuous distribution of the deformation field, effectively improving spatial resolution and data integrity. The generated multi-scale characteristic spectrum can comprehensively reflect the complex behavior of the deformation field, providing a refined multi-level data foundation for subsequent analysis. It ensures physical consistency during the deformation reconstruction process, enhancing the scientific nature and reliability of the results. It achieves high-precision conversion from discrete points to a continuous field, improving the ability to describe the distribution characteristics of the deformation field. It provides an important reference for deformation analysis, engineering monitoring, and structural optimization design of rockfill dams, and has high practical engineering application value.

[0115] Traditional interpolation methods are usually based on pure mathematical models and lack physical constraints. Therefore, according to one aspect of the present application, the steps of performing variational mode decomposition, generating multi-scale characteristic spectra and constructing multi-scale kernel functions, and establishing a set of physical response kernel functions specifically include:

[0116] The variational modal decomposition algorithm is applied to the deformation monitoring data in the standardized monitoring data set to decompose the deformation field signal into modal components with different frequency characteristics, identify multi-level characteristic scales, and generate a multi-scale characteristic spectrum;

[0117] The variational mode decomposition algorithm is: min{uk,ωk}{Σk||∏t[(Δ(t)+j / πt)*uk(t)]·e -jωkt || 2 2}; st Σk uk(t) = u(t); where uk(t) represents the modal component of the kth scale; ωk represents the center frequency of the kth mode; Δ(t) represents the Dirac function; j represents the imaginary unit; * represents the convolution operator; u(t) represents the original signal; ∏ is the partial derivative; t is the time variable; e -jωkt is a complex exponential function; ||·|| 2 2 is the square of the L2 norm; k is the modal index; ∏t is the time partial derivative operator;

[0118] According to the characteristic functions of mechanical differential operators in the set of deformation physical constraints, a set of physical response kernel functions that satisfy the mechanical equilibrium equations is constructed;

[0119] The physical response kernel function set is constructed as follows: K(x, y) = Σn[Ξn(x)·Ξn(y) / λn]; where Ξn(x) represents the nth eigenfunction of the mechanical differential operator L, satisfying L[Ξn(x)] = λn·Ξn(x); λn represents the corresponding eigenvalue; K(x, y) represents the physical response kernel function; x and y represent the spatial coordinates; and n represents the eigenmode index. The combination of variational mode decomposition and the physical response kernel function organically integrates the multiscale characteristics of the signal with physical constraints, improving the physical rationality of the deformation field reconstruction.

[0120] According to one aspect of the present application, the steps of performing multi-scale kernel function interpolation calculation and outputting a high-precision continuous distribution of the deformation field include:

[0121] Apply the physical response kernel function set to the deformation monitoring data in the standardized monitoring data set, establish an interpolation equation system, and output a high-precision continuous distribution of the deformation field;

[0122] Establish an interpolation equation system: u(x) = Σi[αi·K(x,xi)] + p(x); Σi[αi·p(xi)]= 0 forall p∈P; where the coefficient αi must satisfy: Σj[αj·K(xi,xj)] + p(xi) = ui for i=1,2,...,N; u(x) represents the deformation field reconstructed by interpolation; xi represents the position of the i-th monitoring point; ui represents the monitoring value of the i-th monitoring point; K(x,y) represents the physical response kernel function; αi represents the undetermined coefficient; p(x) represents a polynomial function to ensure the completeness of the interpolation; P represents the polynomial space of a given order; N represents the total number of monitoring points; i and j are monitoring point indices; x is an arbitrary spatial point; for different spatial regions, the kernel function parameters are adaptively adjusted to optimize the interpolation accuracy; the cross-validation method is used to evaluate the accuracy of the interpolation results to form a continuous distribution of the deformation field with a confidence index;

[0123] This embodiment ensures that the deformation field satisfies the mechanical equilibrium condition while ensuring the accuracy of data fitting, thereby improving the accuracy and physical rationality of deformation field reconstruction.

[0124] According to one aspect of the present application, the steps of reconstructing discrete monitoring point data by a multi-scale physical constraint kernel function interpolation method to generate a continuous distribution of the deformation field include:

[0125] Read deformation monitoring data and rockfill state field, construct a multi-resolution hierarchical structure, perform spectral decomposition and harmonic residual calculation, and obtain anomaly point set and anomaly credibility index; read anomaly point set, deformation monitoring data and multi-scale spectral residual set, estimate the preliminary repair value of the anomaly point, establish a physical constrained optimization problem, and obtain a constrained optimization problem; read the constrained optimization problem, apply the iterative harmonic solution to solve it, and perform uncertainty assessment and verification to obtain repaired deformation monitoring data and repair uncertainty index; read the repaired deformation monitoring data, apply the physical response kernel function interpolation calculation, and obtain a high-precision continuous distribution of the deformation field.

[0126] The construction of multi-resolution hierarchical structure includes:

[0127] Read deformation monitoring data and standardize them to obtain a standardized measurement matrix and correlation index. Perform quadtree decomposition based on the spatial distribution of monitoring points: initialize the quadtree root node, which contains all monitoring points. If the number of points contained in a node is greater than the threshold and the area size is greater than the minimum size, it is divided into four child nodes. Recursively perform the above steps until the termination condition is met. Define the resolution level according to the quadtree structure: H1 = {the center point of each node in the first layer of the tree}, H2 = {the center point of each node in the second layer of the tree}, ..., H m = {the center point of each node in the deepest level of the tree}, where the center point is defined by the weighted average of the monitoring points in the area, with the weights based on the measurement reliability.

[0128] Among them, spectral decomposition and residual calculation include:

[0129] Read the normalized measurement matrix and the multi-resolution hierarchy, perform singular value decomposition on each resolution level; determine the optimal truncation dimension k_i = argmin_k {k | ∑_{j=1} k σ_j 2 / ∑_{j=1} n σ_j 2 > τ_i}, where τ_i = τ_base + Δτ·(i-1) / (m-1) - ε·PC_i, τ_base is the base threshold, PC_i is the physical constraint satisfaction score of the i-th layer; construct the reconstruction matrix and calculate the spectral residual U_i k = V_i(:,1:k_i) Σ_i(1:k_i,1:k_i) W_i(:,1:k_i) T , SR_i = |U_i - U_i k |. / σ_E_i; calculate the comprehensive anomaly score AS(p) =∑_i w_i SR_i(p) (1 + λ TS(p)), where TS(p) is the temporal stability index and w_i is the hierarchical weight; apply the adaptive threshold to identify outliers τ_AS = μ_AS + η σ_AS (1 + γ D(p)), where D(p) is the isolation index of point p.

[0130] The construction and solution of physical constraint optimization problems include:

[0131] Read the outlier point set and the preliminary repair value set, and construct a physical constraint optimization problem: min_{u*} ||u* - u''|| 2 + μ||L[u*] - f|| 2 + ρ·∑_{q∉O}(u*(q)-u(q)) 2, where u* is the final repair value vector, L is the mechanical differential operator, and f is the volume force term; the adaptive regularization parameter μ(p)=μ0·exp(-λ·CI(p)) is designed, where CI(p) is the anomaly credibility index; the variational multi-grid acceleration strategy is applied to solve: construct a nested space V1 ∽ V2 ∽ ... ∽ V n ; where ∽ represents the true inclusion relation of the set; solve the initial solution on the coarsest grid V1; progressively refine, at each layer V k Initial values are obtained and improved through the previous layer solution and interpolation operators; finally, the repaired deformation monitoring data are obtained; uncertainty assessment is performed, including calculating the residual with physical constraints: r = ||L[u*] - f||; evaluating the spectral residual SR* after repair; and constructing the uncertainty index: U(p) = ω1·r(p) + ω2·SR*(p) + ω3·CI(p).

[0132] According to one aspect of the present application, step S31 is specifically as follows:

[0133] Step S311: Read the raw deformation monitoring data D_raw, which contains the coordinates of the monitoring points {xi, yi, zi} and the corresponding measurement values ui. Perform data format verification, check coordinate consistency, mark significant out-of-limit values, and generate preliminary cleaned data D_clean and a data integrity report R_integrity.

[0134] Read the concrete face rockfill dam completion diagram G and the rockfill state field CSI, and extract the elastic differential operator L = μ▽ 2 (·) + (λ+μ)▽(▽·(·)); where λ and μ are Lame constants, determined based on the material partitioning information. Construct the physical constraints: L[u(x)] = f(x, CSI(x)); B[u(x)] = g(x) on dΩ; where f(x, CSI(x)) is the body force term, related to the permeability-mechanical coupling state indicator CSI(x); B is the boundary operator; g(x) is the boundary condition; and dΩ is the computational domain boundary. Output the deformation physical constraint set C_phys.

[0135] S312. Read the preliminary cleaned data D_clean and construct a deformation measurement matrix U: rows represent time points, and columns represent monitoring points. Normalization is performed using U_norm = (U - μ) / σ, where μ and σ are the mean and standard deviation vectors, respectively. Calculate the temporal correlation R_time and spatial correlation R_space of the measurement matrix. Output the normalized measurement matrix U_norm and the correlation indices {R_time, R_space}.

[0136] Read the preliminary cleaned data D_clean and construct a multi-resolution hierarchy H = {H1, H2,..., H m First, perform quadtree decomposition: initialize the quadtree root node, which contains all monitoring points; if the number of points contained in the node is greater than the threshold T_split and the region size is greater than the minimum size S_min, then divide it into four child nodes; recursively execute until the termination condition is met; then define the resolution level according to the tree structure: H1 = {the center point of each node in the first layer of the tree}; H2 = {the center point of each node in the second layer of the tree}; H m = {the center point of each node in the deepest layer of the tree}; where the center point is defined by the weighted average of all monitoring points in the region, with the weights based on the measurement reliability. Output is the multi-resolution hierarchy H.

[0137] Read the standardized measurement matrix U_norm and the multi-resolution hierarchy H, for each resolution level H i Perform the following steps: Extract the corresponding measurement sub-matrix U_i; Perform singular value decomposition: U_i = V_i Σ_i W_i T ; Determine the optimal truncation dimension k_i = argmin_k {k | ∑_{j=1} k σ_j 2 / ∑_{j=1} n σ_j 2 > τ_i}; where σ_j is the singular value, τ_i is the energy retention threshold, which increases as i increases; construct the reconstruction matrix U_i k = V_i(:,1:k_i) Σ_i(1:k_i,1:k_i) W_i(:,1:k_i) T ; Calculate the reconstruction error E_i = U_i - U_i k Calculate the spectral residual SR_i = |E_i| . / σ_E_i; τ_i = τ_base + Δτ·(i-1) / (m-1) - ε·PC_i; τ_base = 0.85 is the base threshold, Δτ = 0.1 is the layer interval, PC_i is the physical constraint satisfaction score of the i-th layer, and ε = 0.05 is the physical compensation coefficient. Output the multi-scale spectral residual set {SR_i}.

[0138] The multiscale spectral residual set {SR_i} is read, and for each monitoring point p, a comprehensive anomaly score (AS(p) = ∑_i w_i · SR_i(p)) is calculated. w_i is the hierarchical weight, which decreases with increasing resolution. A recursive stability analysis (AS'(p) = AS(p) · (1 + λ · TS(p)) is introduced. TS(p) is the temporal stability indicator, TS(p) = std(AS_t(p)) / mean(AS_t(p)). AS_t(p) is the anomaly score at time t in the past. Outlier identification uses an adaptive threshold: τ_AS = μ_AS + η·σ_AS·(1 + γ·D(p)). D(p) is the isolation index of point p, representing the average distance from neighboring points, with η = 2.5 and γ = 0.3. The final output is the outlier point set O and the anomaly confidence index (CI).

[0139] S313, read the abnormal point set O, the normalized measurement matrix U_norm and the multi-scale spectral residual set {SR_i}, and for each abnormal point p∈O, estimate the preliminary repair value u'(p) = ∑_i v_i·u_i based on the multi-scale reconstruction matrix k (p); where u_i k (p) is the value at point p in the reconstruction matrix of layer i. Vi is the confidence weight, which is inversely proportional to the spectral residual SRi: vi = exp(-α·SRi(p)) / ∑_j exp(-α·SRi(p)). α = 2.0 is the sensitivity parameter. To enhance temporal consistency, a time window smoothing function u''(p) = β·u'(p) + (1-β)·∑_t w_t·u'_t(p) is introduced. Where u'_t(p) is the preliminary repair value at time t in the past, w_t is the time weight, and β = 0.6 is the weight at the current time. The output is the preliminary repair value set u''.

[0140] Read the abnormal point set O, the preliminary repair value set u'' and the deformation physical constraint condition set C_phys, and construct the physical constraint optimization problem: min_{u*} ||u* - u''|| 2 + μ||L[u*] - f|| 2 ; st u*(q) = u(q) for all q belongs to O; where u* is the final repair value vector, u represents the original measurement value, L is the mechanical differential operator, and f is the body force term. To improve computational efficiency, the penalty function method is used to convert it into an unconstrained optimization problem: min_{u*} ||u* - u''|| 2 + μ||L[u*] - f|| 2 + ρ·∑_{q∉O}(u*(q) - u(q)) 2Adaptive regularization parameter μ(p) = μ0·exp(-λ·CI(p)); CI(p) is the anomaly confidence index, μ0=0.5 is the base parameter, and λ=2.0 is the attenuation coefficient. This allows the algorithm to rely more on physical constraints than preliminary repair values for highly anomalous points. Output is the constrained optimization problem P_opt.

[0141] Read the constrained optimization problem P_opt and apply the iterative harmonic solution to solve it. First, discretize the problem and construct a linear system: (I + μA T A + ρB T B)u* = (u'' + μA T f + ρB T u); where A is the discrete matrix of the mechanical operator L, and B is the selection matrix (to extract non-outlier points). The conjugate gradient method is used to solve the problem, and a variational multi-grid acceleration strategy is introduced: construct a nested space V1 contained in V2 contained in... contained in V n ; Solve the initial solution on the coarsest grid V1; progressively refine, at each layer V k Above, through the front layer solution u*_{k-1} and interpolation operator I k _{k-1} obtains the initial value; performs N conjugate gradient iterations to improve the solution; obtains the final solution u*, and outputs the repaired deformation monitoring data D_fixed.

[0142] The repaired deformation monitoring data D_fixed and the set of physical constraints C_phys are read and uncertainty assessment is performed: the residual with respect to the physical constraints is calculated: r = ||L[u*] - f||; the spectral residual SR* after repair is evaluated; and an uncertainty index is constructed: U(p) = ω1·r(p) + ω2·SR*(p) + ω3·CI(p), where ω1=0.4, ω2=0.3, and ω3=0.3 are weighting coefficients. Cross-validation is performed to assess the repair quality: 50% of non-outlier points are randomly selected as pseudo-outliers; the repair algorithm is applied; the recovery error is calculated; statistical indicators are obtained by repeating the process 20 times; and the repair uncertainty index U and the repair quality assessment report R_quality are output. Physical constraints are introduced in the anomaly detection phase to assist in spectral truncation assessment. In the repair phase, preliminary statistical repair and physical constraint optimization are unified into a single framework.

[0143] According to one aspect of the present application, the steps of generating a comprehensive report on fracture properties through fracture feature extraction and visual decision support include:

[0144] S41, reading the panel stress field and panel deformation field, calculating the crack sensitivity index of each area of the panel based on crack mechanics theory, and outputting a crack sensitive area map;

[0145] S42, reading underwater detection data and a crack sensitive area map, applying image enhancement technology and feature extraction algorithms tailored to underwater environments, identifying crack locations, directions, and widths, and outputting crack feature data;

[0146] S43, reading crack characteristic data, panel stress field, and panel deformation field, comprehensively analyzing the cause, development status, and risk level of the cracks, and outputting the crack property determination result;

[0147] S44. Read the crack property determination results, crack characteristic data and crack sensitive area map, generate a multi-dimensional crack property visualization expression, provide targeted maintenance suggestions based on risk analysis, and output a comprehensive crack property report.

[0148] According to one aspect of the present application, physical trigger condition association analysis and adaptive threshold transition point identification include:

[0149] Read the physical parameters in the standardized monitoring data set, extract the effective stress and pore water pressure, and calculate the rate of change of physical parameters: Rσ(x,t) = |σ'(x,t) - σ'(x,t-Δt)| / σ'_max, Rp(x,t) = |p(x,t) - p(x,t-Δt)| / p_max; construct the physical trigger correlation function C_j(x,t) = max(corr(R'_j(x,t±τ), Rσ(x,t)), corr(R'_j(x,t±τ), Rp(x,t))), where corr represents the sliding window correlation coefficient and τ is the time lag / lead range; design the adaptive threshold function τ_j(x,t) = μ_j(x) + α·σ_j(x,t)·(1 + β·exp(-γ·t_s(x)))·(1 - ω·C_j(x,t)), where μ_j(x) is the base threshold, σ_j(x,t) is the local standard deviation, and t_s(x) is the time since the last transition point; identify the set of transition points and perform multi-scale fusion TP_j(x) = {t | R'_j(x,t)>τ_j(x,t) and local maximum}, TP(x) = Filter({TP_j(x)}, d_min, n_min).

[0150] Among them, multi-scale hierarchical physical modeling includes:

[0151] Construct a low-frequency physical diffusion model: da_J(x,t) / dt = ▽·[D_L(x,φ(x,t))·▽a_J(x,t)]+ S_L(φ(x,t)), where D_L(x,φ) = D_0·exp(-α_e·σ'(x,t))·exp(λ·a_J(x,t))·R(θ(x)), S_L(φ) = r_L·(a_J,eq(φ) - a_J); construct a detail scale adaptive model: dc_j(x,t) / dt =D_j·▽ 2 c_j(x,t) + v_j·▽c_j(x,t) + S_j(φ,a_J,t), where S_j(φ,a_J,t) = S_j,eq(φ,a_J) + S_j,trans(t) + S_j,noise, S_j,trans(t) = ∑_{t_c∈TP} A_j(t_c)·h_j(t-t_c), h_j(t) = (1-exp(-ρ_j·t)) / (1+exp(-ρ_j·t)); establish a cross-scale coupling mechanism: F_j(t) = ∑_i I_{i,j}·∫ K_{i,j}(t-τ)·c_i(τ)dτ, dc_j(x,t) / dt = ... + F_j(t), and at the same time modify a_J,eq(φ) to a_J,eq(φ,{c_j}) to form a bidirectional coupling.

[0152] Among them, consistency optimization based on physical constraints and Bayesian posterior analysis include:

[0153] Construct the physical consistency constraint optimization problem min_{a',c'} ||a'-a|| 2 + ∑_j||c'_j-c_j|| 2 + λ||L[CSI'] - f(φ)|| 2 , where λ(x) = λ_0·(1 + κ·Ξ(x)), Ξ(x) is a physical importance indicator; Bayesian posterior analysis is applied p(θ|D) ∝ p(D|θ)p(θ), where θ is the model parameter vector and D is the historical observation data; the Markov chain Monte Carlo method is used to sample the posterior distribution and generate the parameter sample set {θ i}; Perform model prediction for each sample and calculate the predicted mean and variance: CSI_mean(x,t) = (1 / N)·∑_i CSI i (x,t), CSI_var(x,t) = (1 / N)·∑_i (CSI i (x,t) - CSI_mean(x,t)) 2; Calculate the prediction reliability index RI(x,t) = (1 - U_pred(x,t) / CSI_pred(x,t))·(1 - RMSE / CSI_range).

[0154] According to one aspect of the present application, step S42 may be:

[0155] S421. Read the original image I of the underwater crack. For each pixel (x, y), calculate the mean μW(x, y) and standard deviation σW(x, y) of a local window W(x, y) centered at that pixel. The window size is 15 × 15 pixels. Apply the local illumination normalization formula I_norm(x, y) = (I(x, y) - μW(x, y)) / σW(x, y) to obtain the illumination normalized image I_norm.

[0156] Read the illumination-normalized image I_norm and separate its three RGB channels: IR, IG, and IB. Apply channel-wise analysis to calculate the contrast scores (CR, CG, and CB) and signal-to-noise ratios (SR, SG, and SB). Construct a crack enhancement weight vector w = [wR, wG, wB]: wR = 0.2 + 0.4 CR + 0.4 SR; wG = 0.3 + 0.4 CG + 0.3 SG; wB = 0.5 + 0.3 CB + 0.2 SB. Normalize the weight vector and perform weighted fusion: I_fused = wR IR + wG IG + wB IB, yielding the fused grayscale image I_fused.

[0157] S422. Read the fused grayscale image I_fused and apply a four-layer discrete wavelet transform, selecting the sym4 wavelet basis function. Calculate the wavelet coefficients at four scales: the low-frequency approximation coefficient cA_4 and the detail coefficients {cH_j, cV_j, cD_j}, where j = 1, 2, 3, and 4, representing horizontal, vertical, and diagonal details, respectively.

[0158] Read the horizontal detail coefficients cH_j and vertical detail coefficients cV_j, and for each scale j and spatial position (x,y), calculate the directional coherence measure M_j: M_j(x,y) = |∑ω cH_j(x,y)·cH_j(x+ω,y) + cV_j(x,y)·cV_j(x,y+ω)| / sqrt(cH_j(x,y) 2 + cV_j(x,y) 2); where ω∈[-3,3] represents the local neighborhood range. To improve sensitivity to small cracks, a directional consistency compensation γ_j(x,y) = max(0, 1 - var(θ(x+Δx,y+Δy)) / π) is introduced; where θ represents the local gradient direction and var represents the variance. The resulting compensated directional coherence measure M'_j(x,y) = M_j(x,y)·(1 + α_j·γ_j(x,y)); where α_j is the scale adaptation coefficient, α_1=0.8, α_2=0.6, α_3=0.4, and α_4=0.2. The output directional coherence measure set {M'_j} is obtained.

[0159] Read the detail coefficients {cH_j, cV_j, cD_j} of each scale and convert them into amplitude and phase representation: A_j(x,y) =sqrt(cH_j(x,y) 2 + cV_j(x,y) 2 + cD_j(x,y) 2 ); Ξ_j(x,y) = atan2(cV_j(x,y), cH_j(x,y)); Cross-scale phase consistency is calculated: PC(x,y) = |∑j e_j·A_j(x,y)·cos(Ξ_j(x,y) – Ξ*(x,y))| / ∑j A_j(x,y); where e_j is the scale weight, e_1=0.4, e_2=0.3, e_3=0.2, e_4=0.1, and Ξ*(x,y) is the weighted average phase. To improve adaptability to underwater environments, a brightness compensation function η(x,y) = 1 / (1 +exp(-k·(cA_4(x,y) - cA_th))) is introduced; where k=5 and cA_th is the adaptive threshold, determined from cA_4 using the Otsu method. Finally, the adjusted phase consistency is obtained: PC'(x,y) = PC(x,y)·η(x,y), and the phase consistency diagram PC' is output.

[0160] S423. Read the directional coherence measure set {M'_j} and the phase consistency map PC', and construct the adaptive enhancement function f_j(M',PC') = (1 + tanh(α_j·(M'_j - β_j))) / 2 · (1 + γ_j·PC'); where α_j is the sensitivity parameter, β_j is the threshold parameter, and γ_j is the phase enhancement coefficient. Adaptive determination of the parameters: α_j = α_0·(1 +Δ·SNR_j); β_j = β_0·(1 - ρ·SNR_j); γ_j = γ_0·(1 + ω·PC_mean); where SNR_j is the estimated j-th scale signal-to-noise ratio, PC_mean is the average phase consistency, α_0=3.5, β_0=0.45, γ_0=0.6, Δ=0.3, ρ=0.2, ω=0.5. Output enhancement function set {f_j}. Read wavelet coefficients {cH_j, cV_j, cD_j} and enhancement function set {f_j}, and perform selective enhancement: cH'_j(x,y) = cH_j(x,y)·f_j(M'_j(x,y), PC'(x,y)); cV'_j(x,y) = cV_j(x,y)·f_j(M'_j(x,y), PC'(x,y)); cD'_j(x,y) = cD_j(x,y)·f_j(M'_j(x,y), PC'(x,y));

[0161] To prevent over-enhancement, an amplitude limiting function cH''_j = sign(cH'_j)·min(|cH'_j|, λ_j·std(cH_j)) is introduced; where λ_j is the scale-adaptive coefficient, and its values for j=1, 2, 3, and 4 are 4.0, 3.5, 3.0, and 2.5, respectively. The vertical and diagonal coefficients cV''_j and cD''_j are processed identically. The low-frequency coefficients remain unchanged: cA'_4 = cA_4. The enhanced wavelet coefficients {cA'_4, cH''_j, cV''_j, cD''_j} are output. The enhanced wavelet coefficients {cA'_4, cH''_j, cV''_j, cD''_j} are read and an inverse wavelet transform is performed to reconstruct the image I_rec. To further optimize crack details, local statistical enhancement is applied to the reconstructed image: I_enh(x,y) = I_rec(x,y) + λ·(I_rec(x,y) - μ_L(x,y))·(σ_G / σ_L(x,y)). Where μ_L and σ_L are the local mean and standard deviation, σ_G is the global standard deviation, and λ = 0.7 is the enhancement factor. This outputs the enhanced crack image I_enh.

[0162] S424. Read the enhanced crack image I_enh, and obtain the preliminary crack binary image B_init by applying adaptive threshold segmentation. Construct a directional linear structural element SE_θ, where θ ∈ {0°, 45°, 90°, 135°} and the length is 7 pixels. Perform morphological opening operations for each direction, and then take the minimum response B_open(x,y) = min_θ{(B_init O SE_θ)(x,y)}; where O represents the morphological opening operation. Then apply the conditional dilation operation B_dilate = Δ_I_enh(B_open, T); where Δ_I_enh represents the dilation operation conditional on I_enh, and T is the termination condition. Output the crack binary image B_dilate.

[0163] Read the crack binary image B_dilate and the enhanced crack image I_enh, and obtain the crack centerline S by applying the skeleton extraction algorithm. For each point p on the centerline S, calculate the crack width W(p) perpendicular to the direction θ(p) as W(p) = argmin_d{I_enh(p+d·n(p))<T_edge and I_enh(p-d·n(p))<T_edge}; where n(p) is the normal vector at point p, and T_edge is the edge threshold. Calculate the crack length L as the total length of the centerline S. Integrate to obtain the crack feature dataset F = {S, W, L,θ}, which includes the centerline, width distribution, total length, and direction distribution. Read the crack feature dataset F and the enhanced crack image I_enh, and construct a false color enhancement display. Map the crack width W to the color space C(p) = ColorMap(W(p) / W_max); where ColorMap is a color mapping from blue (thin) to red (thick). Generate the enhanced visualization image I_vis = α·I_enh + (1-α)·Rasterize(S, C); where Rasterize represents rendering the centerline S as a raster image according to the color C, and α = 0.4 is the mixing coefficient. Output the crack visualization image I_vis.

[0164] In this embodiment, a parameter adaptive adjustment method based on local signal-to-noise ratio and phase characteristics is established, enabling the system to adapt to the imaging condition differences in different underwater regions. Nonlinear selective enhancement based on crack morphological features is implemented in the wavelet domain, solving the traditional contradiction between "over-smoothing" and "noise retention".

[0165] In a specific embodiment of the present application, a method for identifying underwater crack properties of a panel rockfill dam is suitable for the safety monitoring and maintenance of large panel rockfill dams. The overall process includes four main steps: multi-source data acquisition and preprocessing, rockfill body permeability-mechanical coupling state modeling, panel deformation and stress analysis, and crack feature extraction and visualization decision support. The implementation process of this method will be described in detail below, combined with actual monitoring data from a large panel rockfill dam. It should be noted that non-innovative points such as knowledge or skills known to those skilled in the art are briefly described.

[0166] Step 1: Multi-source data collection and preprocessing.

[0167] 1.1. Acquisition of multi-source monitoring data.

[0168] In this example, the original monitoring data of the concrete face rockfill dam were first obtained, including: 3D geometric data of the above-water portion of the concrete face, the dam crest, and the downstream dam surface obtained by 3D laser scanning, with a spatial resolution of 5 mm; underwater sonar scanning data, including underwater dimensions and blanket morphology information of the concrete face, with a resolution of 10 mm; dam body deformation data obtained by a water-tube sedimentation meter, with a measurement accuracy of 0.1 mm; stress and strain monitoring data, including stress and strain time series data from 25 monitoring points; and pore water pressure monitoring data from 32 monitoring points inside the dam body.

[0169] 1.2. Data quality assessment and anomaly detection.

[0170] Statistical methods were used to identify and flag outliers: a 3σ standard was calculated for each monitoring point, and data outside the mean ±3σ range were marked as potential outliers. Further verification was performed using the Grubbs test, calculating the statistic G = |outlier - mean| / standard deviation. A data quality index (DQI) was calculated for each monitoring point = (1 - proportion of outliers) × (1 - proportion of missing values) × (1 - standard deviation / mean). This step identified a total of 57 outliers and 12 monitoring points with data quality below 0.75.

[0171] 1.3. Spatiotemporal alignment and fusion preprocessing.

[0172] To address the spatiotemporal inconsistency issues of multi-source data: a unified coordinate system is established based on high-precision GPS coordinates and the measurement control network; the Helmert seven-parameter transformation is applied to convert data from different sources into a unified coordinate system; data time is aligned based on timestamps, and linear interpolation is used to handle different sampling frequencies; a spatiotemporal aligned dataset is generated to achieve consistent expression of multi-source data.

[0173] 1.4. Adaptive noise reduction and data enhancement.

[0174] According to different sensor characteristics and environmental interference: wavelet denoising algorithm is applied to process noise data, and db4 wavelet basis function is selected; for stress monitoring data, threshold coefficient λ = 3.5×sqrt(2logN), N is the data length; for underwater sonar data, adaptive Bayesian shrinkage threshold method is applied, and threshold calculation is λj = σ 2 / σj,σ 2 is the noise variance, and σj is the signal variance; the comprehensive data reliability index is calculated: RI = w1×DQI + w2×SNR + w3×TC, where DQI is the data quality index, SNR is the signal-to-noise ratio, and TC is the temporal consistency, with weights w1=0.5, w2=0.3, and w3=0.2 respectively; through this step, the noise level is reduced by 85% and the data reliability is improved by 32%, laying the foundation for subsequent analysis.

[0175] This embodiment combines spatiotemporal alignment of multi-source data with adaptive noise reduction, addressing the difficulty of unified processing of heterogeneous multi-source data in traditional methods. In particular, the adaptive wavelet noise reduction algorithm dynamically adjusts the threshold based on the characteristics of different sensors, improving noise reduction efficiency by approximately 30% compared to traditional fixed-threshold methods while effectively preserving important feature information. The introduction of a comprehensive data reliability index provides a data quality reference for subsequent analysis, enabling a quantifiable assessment of the credibility of analysis results. This is a significant breakthrough in the field of underwater structure inspection.

[0176] Step 2: Modeling the permeability-mechanical coupling state of the rockfill body.

[0177] 2.1. Construction of rockfill physical model.

[0178] Based on the as-built drawings and monitoring data of the concrete face rockfill dam, a 3D geometric model was established, including the structural layers of the rockfill body. A mesh model was generated, with element size varying with depth, ranging from 5 m in the surface area to 10 m in the deep area. Material zoning was determined based on the as-built data, including the transition material area, the main rockfill area, and the cushion layer area. The rockfill body geometry model (containing 10,245 nodes and 45,368 elements) and material zoning information were output.

[0179] 2.2. Construction of permeability-mechanical coupling state index.

[0180] Construct the permeability-mechanical coupling state index (CSI): extract effective stress σ'(x, t), pore water pressure P(x, t) and strain ε(x, t) data from the standardized monitoring data set; calculate the permeability-mechanical coupling state index: CSI(x, t) = α·σ'(x, t) / σ'max + β·P(x, t) / Pmax + γ·ε(x, t) / εmax; where α=0.45, β=0.35, and γ=0.20 are weight coefficients, indicating the relative importance of each factor; σ'max=2.5MPa is the maximum effective stress; Pmax=0.8MPa is the maximum pore water pressure; εmax=0.004 is the maximum strain; x is the spatial coordinate point; and t is the time point.

[0181] Traditional methods typically consider either the stress state or the permeability state independently, making it difficult to characterize the interaction between the two. This example proposes a coupled permeability-mechanical state index (CSI), which unifies three key parameters: effective stress, pore water pressure, and strain, into a single quantitative framework. By assigning different weighting coefficients to reflect the contribution of each factor to the coupled state, this approach overcomes the inability of traditional methods to comprehensively characterize the complex internal state of the rockfill. Practical applications have demonstrated that this index accurately reflects the comprehensive state of the rockfill, improving the accuracy of identifying wetting and stress concentration zones by 42%, providing more reliable baseline data for subsequent analysis.

[0182] 2.3. Stress state sensitive self-organizing map.

[0183] Based on the permeability-mechanical coupling state index, a stress-state-sensitive self-organizing map (SOM) was implemented: the input data were the permeability-mechanical coupling state index (CSI) and material partition information; the SOM network parameters were set as follows: a 10×10 neuron grid, an initial learning rate η0=0.1, a time constant τ1=1000, an initial neighborhood radius σ0=5.0, and a time constant τ2=500; the core update formula of the improved SOM algorithm was applied: wij(t+1) = wij(t) + η(t)·h(i, BMU(x), t)·[x - wij(t)]·ψ(σ'(x, t)); where wij(t) is the j-th dimension weight of the i-th neuron at time t; η(t) = η0·exp(-t / τ1) is the learning rate; and h(i, BMU(x), t) = exp(-||ri - rBMU(x)|| 2 / (2σ 2 (t))) is the neighborhood function; σ(t) = σ0·exp(-t / τ2) is the neighborhood radius; ψ(σ'(x, t)) = [1 + exp(-k·(σ'-σ'0))] -1is the stress state sensitivity function; σ'0 = 1.2 MPa is the stress critical value; k = 5.0 is the stress sensitivity parameter. After 5000 iterations, the rockfill state characteristic partition is generated.

[0184] Traditional SOM algorithms have limitations when processing data with nonlinear distribution characteristics. This implementation introduces a stress state sensitivity function, ψ(σ'), which enables the algorithm to differentiate its response sensitivity to different stress state regions, particularly enhancing the identification of critical stress regions. Experimental verification shows that compared with the traditional SOM, the improved algorithm improves classification accuracy by 38% in the critical stress region (0.9 MPa-1.5 MPa) and boundary recognition accuracy by 25%. This improvement is crucial for accurately identifying regions with sudden stress changes within rockfill bodies and provides a more precise spatial zoning basis for subsequent fracture risk assessment.

[0185] 2.4. Dynamic modeling of regional mechanical properties.

[0186] Based on the state characteristics zoning of the rockfill body, a regional mechanical characteristics dynamic model considering stress-seepage coupling is established: the dynamic differential equation is constructed: ∏CSI(x, t) / ∏t = ▽·[D(x, σ', CSI)·▽CSI(x, t)] + Q(x, t, σ'); where D(x, σ', CSI) = D0·exp(-α e σ')·exp(λ·CSI)·R(θ) is the diffusion coefficient tensor; D0=1.5×10 -6 m 2 / s is the reference diffusion coefficient; α e =0.8 MPa -1 is the effective stress influence coefficient; λ = 2.5 is the coupling state sensitivity coefficient; R(θ) is the anisotropic rotation matrix, and the angle between the principal axis and the dam axis is θ = 15°; Q(x, t, σ') is the source-sink term, which represents the state change rate caused by local stress change, Q(x, t, σ') = β1·∏σ' / ∏t·[1-exp(-β2·σ')], where β1 = 0.01 and β2 = 1.2 MPa -1 The finite element method was used to solve the dynamic equations with a time step of 1 day to generate a set of local mechanical response models.

[0187] Traditional seepage-mechanical analysis typically uses a fixed diffusion coefficient, which cannot reflect the complex nonlinear characteristics of the rockfill. This example expresses the diffusion coefficient as a function of the effective stress and the coupled state index (D(x, σ', CSI), and introduces a source-sink term (Q(x, t, σ')) to characterize local state changes. This representation accurately describes the nonlinear mechanical response of the rockfill under complex stress states, particularly capturing the impact of stress changes on permeability. Field data verification demonstrates that compared to traditional fixed-coefficient models, this dynamic model improves the prediction accuracy of internal state changes in the dam by 45% and the sensitivity of identifying local anomaly areas by 57%, providing a reliable basis for accurately assessing crack risk areas.

[0188] 2.5. Boundary adaptive fusion.

[0189] Adopt boundary adaptive fusion method to achieve seamless connection of different partition models: boundary adaptive fusion algorithm: CSI(x) = Σ i [W i (x)·CSI i (x)] W i (x) = R i ·ρ i (x) / Σ j [R j ·ρ j (x)] ρ i (x) = exp(-d i (x) 2 / 2l i 2 ); where CSI(x) is the fused state indicator field; CSI i (x) is the state index of the i-th local model; W i (x) is the weight of the i-th model at position x; R i is the reliability score of the i-th model, ranging from [0.75-0.95], determined based on the data reliability index; ρ i (x) is the distance decay function; d i (x) is the distance from point x to the center of the ith region; l i is the characteristic length scale of the i-th region, which is determined adaptively according to the partition size and ranges from [10m-30m]. i , achieving smooth transition of the boundary area; outputting the spatial distribution characteristics of the rockfill state, including the global CSI field and regional boundary transition characteristics.

[0190] Traditional regional models often have discontinuities at the boundaries, resulting in jumps in the calculation results. The boundary adaptive fusion method proposed in this embodiment achieves seamless connection of different partition models by introducing a weight function based on data reliability and spatial distance. In particular, by adaptively adjusting the feature length scale l i , enabling the model to automatically optimize the width of the boundary transition region based on local characteristics. Practical applications have demonstrated that this embodiment eliminates discontinuities at model boundaries, resulting in smoother and more natural state field transitions. This improves the accuracy of boundary region state prediction by 36%, providing high-quality state space distribution features for subsequent deformation and stress field calculations.

[0191] 2.6. Simulation of dynamic evolution of time-varying states and quantification of uncertainty.

[0192] Based on the spatial distribution characteristics of the rockfill state, a dynamic evolution simulation of the time-varying state was performed: the time-varying differential equations were solved through prediction-correction iterations with a time step of 1 day; historical monitoring data were used to calibrate the model parameters to optimize the accuracy of evolution prediction; the uncertainty of the state field was quantified through Monte Carlo simulation and variational Bayesian inference: 1000 Monte Carlo simulations were performed, and the input parameters were randomly perturbed by ±10%; the standard deviation distribution of the state field was calculated to generate a state uncertainty field; high uncertainty areas were identified to provide a reference for risk assessment; and the rockfill state field (with a resolution of 1m×1m×1m) and the state uncertainty field were output.

[0193] Traditional state evolution analysis often ignores the impact of parameter uncertainty, resulting in an incomplete risk assessment. This embodiment combines Monte Carlo simulation and variational Bayesian inference to not only predict the evolution of the state field but also quantify the uncertainty distribution of the prediction results. This can identify high-risk areas and provide more comprehensive information support for decision makers. Practical applications have shown that the average relative error of state evolution predictions has been reduced to 6.3%, a 42% reduction compared to traditional methods. At the same time, uncertainty quantification has increased the early identification rate of high-risk areas by 65%, enhancing risk warning capabilities.

[0194] Step 3: Panel deformation and stress analysis.

[0195] 3.1. Interpolation and reconstruction of discrete monitoring point data.

[0196] 3.1.1. Physical constraint extraction

[0197] Based on the as-built drawings of the concrete face rockfill dam and the state field of the rockfill body: extract the deformation physical constraint conditions based on the principle of mechanical equilibrium: L[u(x)] = f(x, CSI(x)) B[u(x)] = g(x) on ∏Ω; where L = μ▽ 2(·) + (λ+μ)▽(▽·(·)) is the elastic differential operator; u(x) is the displacement field; f(x, CSI(x)) is the body force, which is related to the permeability-mechanical coupling state indicator CSI(x); B is the boundary operator; g(x) is the boundary condition; λ = 25 GPa and μ = 12 GPa are the Lame constants, determined based on material properties; ∏Ω is the computational domain boundary. Generate a set of deformation physical constraints, including equilibrium equations and boundary conditions.

[0198] 3.1.2. Hierarchical feature scale analysis.

[0199] Apply variational mode decomposition to deformation monitoring data: Variational mode decomposition algorithm: min{uk,ωk}{Σk||∏t[(Δ(t)+j / πt)*uk(t)]·e -jωkt || 2 2}; st Σk uk(t) = u(t); where uk(t) is the modal component of the kth scale; ωk is the center frequency of the kth mode; Δ(t) is the Dirac function; j is the imaginary unit; * is the convolution operator; e -jωkt is a complex exponential function; ||·|| 2 2 is the square of the L2 norm; u(t) is the original signal; t is the time variable; k is the modal index; and ∏t is the time partial derivative operator. The deformation field is decomposed into six main modes, covering deformation characteristics at different time scales. The output is a multi-scale characteristic spectrum containing frequency and energy distribution information.

[0200] Traditional deformation analysis typically treats deformation as a single-scale phenomenon, making it difficult to distinguish deformation characteristics at different time scales. This embodiment introduces variational modal decomposition technology to decompose the deformation field signal into modal components with different frequency characteristics, enabling multi-scale feature extraction. Experimental results show that this embodiment successfully identifies seasonal variations in dam deformation (low-frequency modes), water level changes (medium-frequency modes), and local structural responses (high-frequency modes). Compared with traditional Fourier analysis, feature recognition accuracy is improved by 53%, providing a scientific basis for subsequent kernel function construction.

[0201] 3.1.3. Construction of physical response kernel function.

[0202] Combining the set of deformation physical constraints and the multiscale characteristic spectrum: The physical response kernel function construction method is: K(x, y) = Σn[Ξn(x)·Ξn(y) / λn]; where Ξn(x) is the nth eigenfunction of the mechanical differential operator L, satisfying L[Ξn(x)] = λn·Ξn(x); λn is the corresponding eigenvalue; K(x, y) is the physical response kernel function; x and y are spatial coordinates; n is the characteristic mode index, taking the first 50 main modes. The differential operator characteristic function system is solved using the finite element method; corresponding multiscale physical response kernel functions are constructed for different scale characteristics; and the output physical response kernel function set contains kernel functions at six scales.

[0203] Traditional kernel function interpolation methods are typically based on purely mathematical models (such as radial basis functions) and lack physical constraints. This embodiment constructs a physical response kernel function that satisfies the mechanical equilibrium equations. By solving the characteristic function system of mechanical differential operators, it ensures that the interpolation results conform to physical laws. Practical verification has shown that this physical response kernel function improves the physical rationality of the interpolation results by 78% while maintaining data fitting accuracy. In particular, prediction accuracy increased by 65% in areas with sparse monitoring points, providing theoretical support for the accurate reconstruction of the panel deformation field.

[0204] 3.1.4. Multi-scale kernel function interpolation calculation.

[0205] The physical response kernel function set is applied to process the deformation monitoring data: the multi-scale kernel function interpolation equation system is as follows: u(x) = Σi[αi·K(x,xi)] + p(x); Σi[αi·p(xi)] = 0 for all p∈P; where u(x) is the deformation field reconstructed by interpolation; xi is the position of the i-th monitoring point; K(x,y) is the physical response kernel function; αi is the unknown coefficient, which is determined by solving the equation system: Σj[αj·K(xi,xj)] + p(xi) = ui for i=1,2,...,N; ui is the monitoring value of the i-th monitoring point; p(x) is a second-order polynomial function to ensure the completeness of the interpolation; P is the polynomial space of a given order; N is the total number of monitoring points, which is 168 in this case. The kernel function parameters are adaptively adjusted for different spatial regions. Ten-fold cross-validation is used to evaluate the interpolation accuracy, and the root mean square error (RMSE) and consistency index (CI) are calculated. The deformation field is output as a continuous distribution with a spatial resolution of 0.5 m.

[0206] Traditional interpolation methods often suffer from insufficient local accuracy when processing data from non-uniformly distributed monitoring points. The multi-scale kernel function interpolation method proposed in this embodiment achieves high-precision reconstruction of the deformation field by combining physical response kernel functions of different scales. In particular, by adaptively adjusting the kernel function parameters, the algorithm can automatically optimize the interpolation strategy according to the local data density. Practical applications show that the interpolation accuracy (RMSE) of this embodiment is 0.86mm, which is 47% higher than that of the traditional radial basis function method; in areas with sparse monitoring points, the prediction accuracy is improved by 62%, enhancing the reliability of deformation field reconstruction.

[0207] 3.2. Mechanical simulation and stress analysis of face plate-rockfill contact.

[0208] Based on the reconstructed deformation field and rockfill state field: a panel-rockfill contact mechanical model is constructed, and the augmented Lagrangian method is used to deal with contact constraints; the influence of the permeability-mechanical coupling state on the contact characteristics is considered, and the contact stiffness and friction coefficient are modified; the stress distribution on the contact surface, including normal stress and tangential stress, is calculated; based on the contact stress distribution and panel structure information, the stress state of the panel under complex loads is analyzed; the panel stress field is solved using the finite element method with a mesh size of 0.5m and a total of 12,568 elements; the panel stress field is output, including principal stress, shear stress, and equivalent stress.

[0209] Traditional contact analysis typically uses fixed contact parameters, making it difficult to reflect the impact of rockfill conditions on contact characteristics. This embodiment incorporates a permeability-mechanical coupling state into the contact model, modifying the contact stiffness and friction coefficient to achieve a more accurate contact mechanics simulation. Comparative analysis shows that the contact model incorporating the coupling state improves prediction accuracy by 39% compared to the traditional model. In particular, stress prediction errors in stress concentration areas are reduced by 56%, providing more reliable stress field data for crack risk assessment.

[0210] 3.3. Deformation prediction and verification.

[0211] Combining the rockfill state field, state uncertainty field, and historical deformation data: a deformation prediction model based on a long short-term memory network (LSTM) was constructed; the model inputs included factors such as the permeability-mechanical coupling state indicator (CSI), water level changes, and temperature; the model was trained using a historical dataset containing five years of monitoring data, of which 70% was used for training and 30% for validation; prediction uncertainty was assessed through Monte Carlo random sampling; and a verified panel deformation field and prediction accuracy evaluation report were output, with an average relative error of 4.2%.

[0212] Step 4: Crack feature extraction and visualization decision support.

[0213] 4.1. Identification of crack-sensitive areas.

[0214] Based on the panel stress field and deformation field: Calculate the crack sensitivity index CFI according to the panel stress state: CFI = (σ1 / ft) α1 + (τmax / fs) α2 + (ε / εc) α3 Where σ1 is the maximum principal stress; τmax is the maximum shear stress; ε is the maximum principal strain; ft = 3.5 MPa is the tensile strength of concrete; fs = 2.8 MPa is the shear strength; εc = 0.0002 is the critical strain; α1 = 2.0, α2 = 1.5, and α3 = 1.0 are weighting indices. Fuzzy comprehensive evaluation is used to determine the crack risk level, which is categorized as low, medium, and high risk. A crack sensitivity map is then output to identify potential crack development areas.

[0215] 4.2. Underwater crack image enhancement and feature extraction.

[0216] Process underwater detection data: Apply adaptive histogram equalization and bilateral filtering to enhance underwater image quality; use an improved U-Net deep learning network for crack segmentation; extract crack centerlines and widths through morphological processing; apply Hough transform to identify the main direction of cracks; output crack feature data, including location, direction, length, and width.

[0217] 4.3. Comprehensive determination of crack properties and visual decision support.

[0218] Comprehensive analysis of crack characteristics: Combined with stress field analysis, determine the cause of cracks: stress type, temperature type, or mixed type; evaluate crack development trends based on stress field, deformation field, and crack characteristics; provide targeted maintenance recommendations based on risk assessment; generate a 3D visualization to display the spatial distribution of cracks and risk level; and output a comprehensive report on crack behavior, including cause analysis, risk rating, and maintenance recommendations.

[0219] Traditional crack analysis is typically limited to surface feature extraction and lacks correlation analysis with internal stress states. This embodiment combines crack characteristics with the permeability-mechanical coupling state, stress field, and deformation field to achieve a comprehensive analysis of crack causes and development trends. In particular, the introduced crack sensitivity index (CFI) comprehensively considers factors such as stress, strain, and material strength, improving the accuracy of crack risk assessment. Practical applications have demonstrated that this embodiment increases the early identification rate of cracks by 58% and the accuracy of risk assessment by 46%, providing a scientific basis and decision-making support for dam safety maintenance.

[0220] This embodiment achieves high-precision identification of underwater crack behavior in face-face rockfill dams through the organic integration of multiple steps. Compared to traditional methods, it offers the following advantages: permeability-mechanical coupled state modeling improves the prediction accuracy of the rockfill's internal state by 45%; multi-scale physical constraint kernel function interpolation improves the accuracy of deformation field reconstruction by 47%; and a contact model that considers the coupled state improves stress prediction accuracy by 39%. This invention achieves high-precision modeling of the permeability-mechanical coupled state of the rockfill, providing a reliable theoretical foundation for subsequent face-face deformation and stress analysis and crack feature extraction, ultimately achieving the goal of accurately identifying underwater crack behavior in face-face rockfill dams.

[0221] The preferred embodiments of the present invention are described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the scope of protection of the present invention.

Claims

1. A method for identifying underwater crack properties of a concrete face rockfill dam, characterized in that: The following steps are involved: Obtain original monitoring data and as-built drawings of the concrete face rockfill dam, and preprocess them to obtain standardized monitoring data sets and data reliability indicators; Based on the standardized monitoring data set, data reliability index and the as-built drawings of the concrete face rockfill dam, the rockfill state field and state uncertainty field are obtained through the rockfill body permeability-mechanical coupling state modeling. Using the rockfill state field, state uncertainty field and as-built drawings of the concrete face rockfill dam, the concrete face deformation field and stress field are obtained through concrete face deformation and stress analysis. Generate a comprehensive report on crack properties based on the panel deformation field, panel stress field, and pre-stored underwater inspection data through crack feature extraction and visual decision support; The process of modeling the rockfill body's permeability-mechanical coupling state includes constructing permeability-mechanical coupling state indicators, specifically: Extract effective stress, pore water pressure, and strain state data from standardized monitoring datasets; Establish a mathematical model that comprehensively considers effective stress, pore water pressure and strain state, and calculate the permeability-mechanical coupling state index CSI(x, t); CSI(x,t) = α·σ'(x,t) / σ'max + β·P(x,t) / Pmax + γ·ε(x,t) / εmax; Where CSI(x, t) represents the permeability-mechanical coupling state index of the spatial point x at time t; σ'(x, t) represents the effective stress; P(x, t) represents the pore water pressure; ε(x, t) represents the strain; σ'max represents the maximum effective stress; Pmax represents the maximum pore water pressure; εmax represents the maximum strain; α, β, and γ represent weight coefficients, and α+β+γ=1; x represents the spatial coordinate point; and t represents the time point.

2. The method according to claim 1, characterized in that The steps of obtaining the rockfill state field and state uncertainty field by rockfill permeability-mechanical coupling state modeling include: Based on the as-built drawings of the concrete face rockfill dam and the structural information in the standardized monitoring data set, a physical model of the rockfill body is established to generate material zoning information. Combining standardized monitoring data sets, data reliability indicators, and material partitioning information, the spatial heterogeneity of the rockfill state is characterized and the spatial distribution characteristics of the rockfill state are formed; Based on the spatial distribution characteristics of the rockfill state and the time series data in the standardized monitoring data set, the dynamic evolution of the time-varying state is simulated and a spatiotemporal evolution model of the rockfill state is constructed; Integrate the spatiotemporal evolution model of rockfill state and data reliability indicators, analyze and quantify uncertainty and its propagation law, and output the rockfill state field and state uncertainty field.

3. The method according to claim 2, characterized in that The steps of characterizing the spatial heterogeneity of the rockfill state and forming the spatial distribution characteristics of the rockfill state include: Extract relevant monitoring data from the standardized monitoring dataset and construct a permeability-mechanical coupling state index that considers effective stress, pore water pressure, and strain state; Based on the permeability-mechanical coupling state index and material partition information, stress state sensitive self-organizing mapping is performed to divide the rockfill state characteristic zones; Based on the zoning of rockfill state characteristics, a regional mechanical characteristic dynamic model considering stress-seepage coupling is established to generate a set of local mechanical response models; Based on the local mechanical response model set and data reliability indicators, boundary adaptive fusion is implemented to achieve seamless connection of different partition models and output the spatial distribution characteristics of the rockfill state.

4. The method according to claim 1, wherein The steps of obtaining the panel deformation field and panel stress field through panel deformation and stress analysis include: Read the deformation monitoring data and rockfill state field in the standardized monitoring data set, reconstruct the discrete monitoring point data through the multi-scale physical constraint kernel function interpolation method, and generate a continuous distribution of the deformation field; Based on the rockfill state field, continuous distribution of deformation field and as-built drawings of the concrete face rockfill dam, a slab-rockfill contact mechanics model is constructed to calculate the contact stress distribution. Based on the contact stress distribution, continuous deformation field distribution, and panel structure information in the as-built drawings of the face rockfill dam, the stress state of the panel under complex loads is analyzed and the panel stress field is calculated. The rockfill state field, state uncertainty field and historical deformation field are continuously distributed, a deformation prediction model is established and verified with measured data, and a verified panel deformation field is output.

5. The method according to claim 4, characterized in that The steps to generate a continuous distribution of the deformation field include: Analyze the as-built drawings of the concrete face rockfill dam and the state field of the rockfill body, extract the deformation physical constraints based on the principle of mechanical equilibrium, and form a set of deformation physical constraints; Perform variational mode decomposition on the deformation monitoring data in the standardized monitoring dataset to identify the multi-level characteristic scales of the deformation field and generate a multi-scale characteristic spectrum; Combining the deformation physical constraint condition set and the multi-scale characteristic spectrum, a multi-scale kernel function that satisfies the physical constraint is constructed, and a physical response kernel function set is established; The physical response kernel function set is used to process the deformation monitoring data in the standardized monitoring data set, perform multi-scale kernel function interpolation calculation, and output a high-precision continuous distribution of the deformation field.

6. The method according to claim 1, characterized in that The steps of crack feature extraction include: Perform local illumination normalization and color channel weighted fusion on the underwater detection data to obtain a fused grayscale image; Apply multi-layer wavelet decomposition to the fused grayscale image to obtain low-frequency approximation coefficients and multi-layer detail coefficients; Calculate the directional coherence measure and phase consistency of the detail coefficients to obtain a directional coherence measure set and a phase consistency map, based on which an adaptive enhancement function is constructed to perform selective enhancement on the wavelet coefficients to obtain enhanced wavelet coefficients; The enhanced wavelet coefficients are subjected to inverse wavelet transform and local statistical enhancement to obtain an enhanced crack image, which is then subjected to morphological sharpening and crack feature measurement to obtain crack feature data.

7. The method according to claim 4, characterized in that The steps for generating a continuous distribution of the deformation field can also be: Read deformation monitoring data and rockfill state field, build a multi-resolution hierarchical structure, perform spectral decomposition and harmonic residual calculation, and obtain anomaly point set and multi-scale spectral residual set; Read the abnormal point set, deformation monitoring data and multi-scale spectral residual set, estimate the initial repair value of the abnormal point, establish a physical constraint optimization problem, and obtain the constrained optimization problem; Read the constrained optimization problem, apply iterative harmonic solution to solve it, perform uncertainty assessment and verification, and obtain repair deformation monitoring data; Read and repair deformation monitoring data, apply physical response kernel function interpolation calculation, and obtain high-precision continuous distribution of deformation field.

8. The method according to claim 7, characterized in that The construction of a multi-resolution hierarchy includes: Read deformation monitoring data and perform standardization processing to obtain a standardized measurement matrix; Perform quadtree decomposition based on the spatial distribution of monitoring points in the standardized measurement matrix, including: Initialize the quadtree root node, which contains all monitoring points; If the number of points contained in the quadtree root node is greater than the threshold and the region size is greater than the minimum size, it is divided into four child nodes; Recursively perform quadtree decomposition until the termination condition is met; Define the resolution hierarchy according to the quadtree structure and output the multi-resolution hierarchy H: Hi = {the center point of each node in the i-th level of the tree}, ..., Hm = {the center point of each node in the deepest layer of the tree}, 1≤i≤m; m is the number of levels of quadtree decomposition; The center point is defined by the weighted average of the monitoring points in the area, with the weights based on the measurement reliability.

9. The method according to claim 2, characterized in that The steps of simulating the dynamic evolution of time-varying states and constructing a spatiotemporal evolution model of rockfill states include: Read the spatial distribution characteristics of the rockfill state and the time series data in the standardized monitoring data set, perform wavelet multi-scale time decomposition, calculate energy analysis and physical trigger condition correlation analysis, and generate a multi-scale wavelet coefficient set and a state transition point set; The multi-scale wavelet coefficient set and state transition point set are read, the physical diffusion model of the low-frequency component and the dynamic model of the medium- and high-frequency detail coefficients are constructed, and a cross-scale dynamic coupling mechanism is established to obtain a multi-scale evolution model set, namely the spatiotemporal evolution model of the rockfill state.

Images