Water conservancy project supervision analysis system and method based on digital twinning

Abstract

This invention discloses a water conservancy project monitoring and analysis system and method based on digital twins, belonging to the field of water conservancy project operation monitoring and management technology. The system includes identifying data points with periodic fluctuation characteristics and sparse spatial distribution, extracting weakly disturbed data subsets, establishing response evolution curves, and identifying potential abnormal areas through neighborhood statistics and clustering. For abnormal areas, it inverts the possible types and distribution directions of micro-damage based on historical time series and environmental boundary conditions to obtain an inverted defect feature set, which is then input into a twin model to perform local thermal-permeability-force coupling calculations to obtain stress response evolution curves within a preset future period. Risk scores are calculated based on the curve's growth rate, step characteristics, and frequency mutations. When the score exceeds a set threshold, the risk level and evolution trend information are output. This invention enables early identification and dynamic evolution analysis of minor hidden dangers, improving the intelligence and precision of water conservancy project safety supervision.

Description

Technical Field

[0001] This invention relates to the field of water conservancy project operation monitoring and management technology, specifically to a water conservancy project monitoring and analysis system and method based on digital twins. Background Technology

[0002] Water conservancy projects, such as dams, sluices, and reservoirs, are crucial infrastructure for controlling water flow, storing and regulating water, flood control and disaster reduction, and irrigation water supply. Their safe operation has a significant impact on regional economic development and the safety of residents' lives and property. Traditional water conservancy project supervision mainly relies on on-site inspections, manual patrols, or monitoring methods based on limited sensor data. These methods suffer from poor timeliness, low perception dimensionality, and delayed fault response, making it difficult to meet the requirements for refined understanding and early warning of equipment operating status under complex conditions.

[0003] Especially under special conditions such as seasonal changes and sudden hydrological shifts, water conservancy projects may experience localized anomalies and hidden dangers such as piping, seepage, and dam cracks. Because these hidden dangers are usually accompanied by extremely minor parameter changes in their initial stages, often below the sensing threshold of equipment, and are irregular and transient in spatial distribution, traditional methods cannot achieve early identification. For example, some earthen dams may experience the potential risk of crack propagation in the dam shoulder section after the dry season and at the beginning of the flood season. Such minute hidden dangers are difficult to identify effectively in sensor responses, becoming a "blind spot" that urgently needs to be addressed in the safety supervision of water conservancy projects.

[0004] Digital twin technology, as an advanced simulation method integrating sensor data, historical operating conditions, engineering modeling, and AI analysis, provides full lifecycle modeling, prediction, and feedback capabilities for complex engineering structures. However, current applications of digital twin systems in water conservancy projects are mostly limited to 3D modeling and status playback, lacking mechanisms for identifying and responding to abnormal features. This makes it impossible to accurately model and predict minor but critical safety hazards, resulting in a "perception gap" between monitoring data and actual safety status. Summary of the Invention

[0005] The purpose of this invention is to provide a water conservancy project monitoring and analysis system and method based on digital twins to address the shortcomings of the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a water conservancy project monitoring and analysis method based on digital twins, comprising:

[0007] Collect multi-source monitoring data during the operation cycle of the target water conservancy project, and synchronize and align the data according to spatial location and time sequence to construct a standardized monitoring dataset D;

[0008] Based on the historical change patterns in dataset D, data points with periodic fluctuation characteristics are identified, and combined with their spatial distribution density, a data subset D1 with weak perturbation characteristics is extracted.

[0009] Multi-parameter joint fitting is performed on the data subset D1 to establish the response evolution curve of each data point, and the curve's slope, fluctuation period, delay time, and local extreme value feature parameters are extracted to form a set of hidden danger feature parameters;

[0010] The set of hazard characteristic parameters is mapped to the spatial coordinate system of the digital twin model, local areas with abnormal evolution trends are screened out, and they are marked as potential abnormal areas in the digital twin model.

[0011] Based on the historical time series data corresponding to the potential anomaly region, a local structural change trend curve of the potential anomaly region is constructed, and combined with the current environmental boundary conditions, the possible micro-damage types and their distribution directions are deduced to obtain the inverted defect feature set.

[0012] The inverted defect feature set is input into the digital twin model, local thermal-permeability-mechanical coupling calculation is performed, and the stress response evolution curve of the potential anomaly region in the future preset period is obtained.

[0013] Based on the growth rate, step characteristics, and frequency mutations in the stress response curve, the risk score V of the potential anomaly region is calculated. If V exceeds the set threshold T, the risk level of the potential anomaly region and the corresponding evolution trend information are output.

[0014] Preferably, the multi-source monitoring data includes dam surface temperature, structural strain, pore water pressure, microseismic amplitude and frequency, and crack features in manually inspected images.

[0015] Preferably, identifying data points with periodic fluctuation characteristics includes:

[0016] For each monitoring point in the standardized monitoring dataset, a time series curve is constructed, and its first derivative is calculated using a sliding time window to extract the local rate of change.

[0017] The time series curve is analyzed in the frequency domain using Fast Fourier Transform to identify the dominant frequency component and calculate the periodicity index P, which is the ratio of the dominant frequency amplitude to the average energy.

[0018] Monitoring points where the periodic index P is greater than the set threshold P0 are identified as data points with periodic fluctuation characteristics.

[0019] Preferably, the multi-parameter joint fitting of the data subset D1 includes:

[0020] For each data point in D1, extract its complete time series and use cubic spline interpolation to fill in the missing data;

[0021] The time series was initially fitted using a multinomial regression model, and the order of the fitting function was determined based on the minimum principle of the Bayesian information criterion.

[0022] Based on regression fitting, the sliding window method is used to further optimize the curve shape locally, extract local change trends, and construct a complete response evolution curve.

[0023] The fitted curve is normalized to standardize the response range of all data points to the [0,1] interval.

[0024] Preferably, the step of extracting the curve's slope, fluctuation period, delay time, and local extremum characteristic parameters includes:

[0025] By calculating the first derivative of the response curve over the entire period, its rate of change function is obtained, thereby extracting the maximum slope value and its corresponding time position, which can be used to reflect the rate characteristics of the response intensity rising phase.

[0026] Periodic analysis was performed on the normalized curve, the autocorrelation function method was used to identify the main period, and the period length was recorded as a characteristic of the fluctuation period.

[0027] The time difference between the initial disturbance point and the maximum response point is calculated using the time delay analysis method and extracted as the delay time parameter.

[0028] Identify the location and amplitude of all local maxima and minima in the curve, and statistically analyze the number of extreme points, average spacing, and amplitude differences to form a subset of local extreme value features.

[0029] Preferably, the screening of local regions with abnormal evolutionary trends includes:

[0030] For each hazard feature vector, bind its corresponding spatial location coordinates (x, y, z) during the data acquisition phase to construct a set of feature points with spatial labels;

[0031] Based on the set of feature points, the fixed radius spherical neighborhood method is used to delineate the local neighborhood of each point, and the mean and variance of the risk characteristics within the neighborhood are calculated.

[0032] Construct a risk density function. If the density of high-risk points in a certain neighborhood exceeds the set threshold ρ0, and the standard deviation of the feature values ​​in the region is lower than the threshold σ0, then it is judged as an abnormal region with a consistent trend.

[0033] Clustering algorithms are used to spatially cluster points that meet certain conditions, and multiple local risk regions with evolutionary consistency are extracted.

[0034] Based on the clustering results, the boundary voxel sets of each region are extracted to form clear spatial boundaries of abnormal regions.

[0035] Preferably, the step of constructing a local structural change trend curve for the potential anomaly region based on the corresponding historical time series data within the potential anomaly region includes:

[0036] Extract time series data of all feature points within the potential anomaly area during the historical operating cycle, including strain, pore water pressure, microseismic response amplitude and frequency parameters;

[0037] For each type of physical quantity, the time series is detrended to eliminate long-term steady-state drift. Then, a smooth response curve is constructed using the moving average and local weighted regression methods.

[0038] By extracting the time-varying gradient, rate of change, and key inflection points of each type of physical quantity through curve fitting, a set of multivariable structural response trend curves is constructed.

[0039] Multiple physical field response curves are aligned along the time axis, and an correlation matrix is ​​established to identify the synergistic relationship between changes in different physical quantities.

[0040] Preferably, the step of deducing the possible types of micro-damage and their distribution direction by combining the current environmental boundary conditions includes:

[0041] Obtain the current environmental boundary conditions corresponding to the abnormal area, including air temperature, water level, reservoir water pressure change rate, rainfall, and seasonal geothermal gradient;

[0042] The correlation analysis between the local structural change trend curve and the current boundary conditions was performed, and the influence intensity was quantified by Pearson correlation coefficient and mutual information index.

[0043] By combining multiple typical response patterns from a pre-defined micro-damage mechanism library, pattern matching is performed on the current trend curve;

[0044] The most likely corresponding micro-damage type is selected by the minimum residual matching principle, and its occurrence direction and spatial offset position are recorded.

[0045] Preferably, the execution of local thermal-permeability-mechanical coupling calculation includes:

[0046] Define the spatial boundary range and time scale of the coupled solution region. The boundary range includes the defect point and its surrounding structural elements with an influence radius of not less than 3 times.

[0047] The current temperature, water pressure, and structural load are imported as initial boundary conditions, and time-step driven boundaries are constructed by combining meteorological and operational forecast data for future preset periods.

[0048] Numerical solution of the heat conduction-seepage-stress coupling equations based on the finite element method or finite volume method;

[0049] At each time step, the principal values ​​of the stress tensor in the local defect region are extracted to form a stress response evolution curve, which is used to characterize the structural evolution trend.

[0050] This invention also provides a water conservancy project monitoring and analysis system based on digital twins, comprising:

[0051] Data acquisition module: Collects multi-source monitoring data during the operation cycle of the target water conservancy project, and synchronizes and aligns the data according to spatial location and time order to construct a standardized monitoring dataset D;

[0052] Data subset construction module: Based on the historical change patterns in dataset D, identify data points with periodic fluctuation characteristics, and extract a set of data subsets D1 with weak perturbation characteristics by combining their spatial distribution density;

[0053] Hazard feature extraction module: Perform multi-parameter joint fitting on data subset D1, establish the response evolution curve of each data point, and extract the curve's slope, fluctuation period, delay time, and local extreme value feature parameters to form a hazard feature parameter set;

[0054] Potential anomaly region marking module: Maps the set of hazard characteristic parameters to the spatial coordinate system of the digital twin model, filters out local areas with abnormal evolution trends, and marks them as potential anomaly regions in the digital twin model;

[0055] Inversion Defect Module: Based on the historical time series data corresponding to the potential anomaly region, construct the local structural change trend curve of the potential anomaly region, and inversely deduce the possible micro-damage types and their distribution directions in combination with the current environmental boundary conditions, thus obtaining the inversion defect feature set;

[0056] Coupled calculation module: Input the inverted defect feature set into the digital twin model, perform local thermal-permeability-mechanical coupled calculation, and obtain the stress response evolution curve of the potential anomaly area in the future preset period;

[0057] Risk level classification module: Based on the growth rate, step characteristics and frequency mutation in the stress response curve, calculate the risk score value V of the potential anomaly area. If V exceeds the set threshold T, output the risk level of the potential anomaly area and the corresponding evolution trend information.

[0058] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0059] 1. This invention extracts the growth rate, step characteristics, and frequency disturbance factors of the stress response curve from multiple dimensions and constructs a weighted risk scoring model, thereby achieving quantitative calculation of structural risks in potential abnormal areas. Compared with traditional monitoring methods that rely on a single indicator (such as absolute stress value or strain rate), this invention can identify risk trends in advance when hidden dangers are still in a small evolutionary stage, effectively avoiding the problem of missed judgments due to insignificant responses, thus significantly improving the sensitivity and foresight of water conservancy project safety supervision.

[0060] 2. This invention constructs a risk scoring threshold system based on historical operational cases and combines it with a digital twin model to dynamically output risk level and trend information, forming a closed-loop mechanism of "calculation-judgment-visualization-feedback". This mechanism not only improves the intuitiveness and operability of risk expression, but also provides data-driven decision-making basis for inspection scheduling, emergency response, and structural reinforcement, which can significantly improve the intelligence and precision of water conservancy project operation safety management. Attached Figure Description

[0061] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0062] Figure 1 This is a mind map of the method of the present invention.

[0063] Figure 2 This is a flowchart of the system modules of the present invention. Detailed Implementation

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

[0065] Example 1, please refer to Figure 1 As shown in this embodiment, the water conservancy project supervision and analysis method based on digital twins includes:

[0066] Collect multi-source monitoring data during the operation cycle of the target water conservancy project, and synchronize and align the data according to spatial location and time sequence to construct a standardized monitoring dataset D;

[0067] Based on the historical change patterns in dataset D, data points with periodic fluctuation characteristics are identified, and combined with their spatial distribution density, a data subset D1 with weak perturbation characteristics is extracted.

[0068] Multi-parameter joint fitting is performed on the data subset D1 to establish the response evolution curve of each data point, and the curve's slope, fluctuation period, delay time, and local extreme value feature parameters are extracted to form a set of hidden danger feature parameters;

[0069] The set of hazard characteristic parameters is mapped to the spatial coordinate system of the digital twin model, local areas with abnormal evolution trends are screened out, and they are marked as potential abnormal areas in the digital twin model.

[0070] Based on the historical time series data corresponding to the potential anomaly region, a local structural change trend curve of the potential anomaly region is constructed, and combined with the current environmental boundary conditions, the possible micro-damage types and their distribution directions are deduced to obtain the inverted defect feature set.

[0071] The inverted defect feature set is input into the digital twin model, local thermal-permeability-mechanical coupling calculation is performed, and the stress response evolution curve of the potential anomaly region in the future preset period is obtained.

[0072] Based on the growth rate, step characteristics, and frequency mutations in the stress response curve, the risk score V of the potential anomaly region is calculated. If V exceeds the set threshold T, the risk level of the potential anomaly region and the corresponding evolution trend information are output.

[0073] The target water conservancy project mainly refers to control-type hydraulic structures such as dams, sluices, and reservoirs, which are primarily made of concrete or earth and rock. Their operating cycle can be divided according to season, hydrological year, or structural maintenance cycle. Generally, the operating cycle is 30 to 365 days, depending on the specific project management requirements.

[0074] During this operating cycle, the collected monitoring data includes, but is not limited to, the following categories:

[0075] Dam surface temperature data: The surface temperature of key parts of the dam is obtained by thermocouples or infrared thermometers, with a sampling frequency of once per hour.

[0076] Structural strain data: acquired by strain gauges or fiber optic sensors embedded in the concrete structure, typically sampled once every 10 minutes;

[0077] Pore ​​water pressure data: Pore pressure gauges are installed inside and outside the dam body or around the foundation to collect pore water pressure that reflects the seepage state. The sampling frequency is the same as that of strain data.

[0078] Microseismic response data: The amplitude and frequency characteristics of the dam body under external disturbances (such as water release, blasting, earthquakes, etc.) are recorded using a ground motion monitoring device, with a sampling frequency of more than 10 times per second;

[0079] Manual inspection image data: Images obtained from on-site manual inspections or drone inspections are processed using image processing algorithms to extract information such as the length, width, and direction of cracks on the structural surface, forming a structural crack feature set.

[0080] Due to the significant differences in the sources and sampling intervals of the aforementioned data, unified processing is required to ensure that the data can be used for evolutionary calculations under the same model. Therefore, this invention employs the following steps to standardize the original data:

[0081] Data preprocessing includes standardizing the time format, removing outliers, and interpolation imputation. All sampling timestamps are converted to Coordinated Universal Time (UTC), and missing measurement points are imputed using local linear interpolation. Outliers that are significantly outside the physical range are replaced with the median.

[0082] Time series synchronization and alignment:

[0083] Considering the inconsistent sampling frequencies of different monitoring quantities, a uniform time step Δt is set, for example, 1 hour. All data are then resampled according to this time step.

[0084] For data with a sampling frequency higher than Δt (such as microseismic data), a moving average is used;

[0085] For data with a sampling frequency lower than Δt (such as artificial images), maintain the most recent value or perform linear prediction within the corresponding time window.

[0086] Spatial location normalization: The spatial locations corresponding to all monitoring data are mapped to a unified three-dimensional coordinate system, based on the dam design drawings or BIM model. The locations of all sensor installation points are represented in the form of "x, y, z", and the crack coordinates of the image data are projected onto the surface of the structural model through image matching.

[0087] The processed data is organized into a dataset D with the following structure. Each data record includes: spatial coordinates (x, y, z); timestamp t; temperature T(x,y,z,t); strain ε(x,y,z,t); pore pressure P(x,y,z,t); microseismic intensity A(x,y,z,t); frequency f(x,y,z,t); and crack feature vector C(x,y,z,t), which includes dimensions such as crack length, width, and orientation.

[0088] The dataset D has a fixed sampling interval in time, and all data indicators have a unified physical unit and dimension.

[0089] In practice, dataset D is stored in a database or multidimensional array format, supporting indexing and batch reading, and is used for subsequent hazard identification, local analysis, and model evolution simulation.

[0090] For various continuous time-series data (such as temperature, pore pressure, strain, etc.) in the monitoring dataset, this invention first performs periodic analysis to identify data points with stable periodic fluctuation behavior. This stage is achieved through the following steps:

[0091] For each sensor or sampling point in the monitoring dataset D, its complete time series within the operating cycle is extracted and denoted as S(t), where t represents time. The local rate of change sequence ΔS(t) is obtained by calculating the difference between adjacent time points, i.e., the first-order difference. This is used to eliminate the influence of the trend term and highlight short-periodic behavior.

[0092] The ΔS(t) sequence is subjected to a Fast Fourier Transform (FFT) to obtain its spectral information. Let the amplitude of the dominant frequency component be A_max, and the total spectral energy be E_total. Then, the periodicity index P at this point is defined as A_max divided by E_total, i.e.:

[0093] Periodicity index P = main frequency amplitude A_max ÷ total frequency domain energy E_total.

[0094] E_total is calculated by the sum of squared energies of all frequency components in the spectrum.

[0095] Set a periodicity identification threshold P0, and select data points with a periodicity index P greater than P0 as "periodic fluctuation feature points". The value of P0 is generally between 0.25 and 0.45, and is specifically determined based on the spectral distribution statistics of historical data.

[0096] For data points selected using periodic indicators, further checks are conducted to determine if there are any missing segments in their time series exceeding a threshold. Only data points with a temporal integrity higher than 95% are retained for subsequent analysis. Ultimately, this stage outputs a set of "periodic fluctuation feature points" for subsequent spatial density analysis.

[0097] Periodic feature points may be concentrated in structurally active regions. To screen out micro-perturbation point sets with spatially sparse distribution characteristics, further three-dimensional spatial analysis is required. The steps are as follows:

[0098] Using the dam's design coordinate system as a reference, the entire monitoring area is divided into a regular cubic grid, with the grid side length set to L, and a recommended range of 2 to 5 meters. Each periodic feature point is assigned to a corresponding grid cell based on its spatial coordinates (x, y, z).

[0099] The number of feature points contained in each grid is counted and denoted as . Calculate the average density N_avg of all grid cells and record the standard deviation σ_N.

[0100] Set the spatial sparsity threshold to N_threshold = N_avg − σ_N. Only retain those that satisfy... Points in the grid with values ​​less than N_threshold constitute a sparse candidate set.

[0101] For each point in the sparse candidate set, calculate its root mean square value of ΔS(t), denoted as RMS_i. Set a perturbation amplitude threshold A0, and retain only data points whose RMS_i is less than A0, forming the final weakly perturbed data subset D1.

[0102] This process ensures that the data points in D1 have "dual sparsity": weak temporal variation and sparse spatial distribution, thus revealing long-term latent risks more representatively.

[0103] To avoid artificially setting overly coarse amplitude limits, this invention constructs an adaptive perturbation amplitude threshold A0 based on statistical methods and data distribution characteristics. The construction method is as follows:

[0104] For the ΔS(t) sequence of all periodic feature points, calculate the standard deviation σ_RMS of its rate of change as a reference for the overall disturbance level.

[0105] The perturbation threshold A0 is defined as k times σ_RMS, i.e., A0 = k × σ_RMS. Here, k is an adjustment factor, and a value between 0.3 and 0.7 is recommended. The specific value should be adjusted according to the application scenario; the lower the k value, the smaller the selected perturbation.

[0106] For the one-dimensional feature set consisting of all RMS values, after calculating its mean and covariance, the Mahalanobis distance formula is used to identify abnormal disturbance points. Points with a Mahalanobis distance greater than a set threshold M0 (e.g., 2.5) will be removed to prevent high-amplitude disturbance points from mistakenly entering the D1 set.

[0107] Verify whether the perturbation distribution of all points in D1 follows a normal distribution. If not, perform a Box-Cox transformation on the RMS value to ensure that the perturbation scale is uniform within a controllable range.

[0108] The disturbance threshold A0 obtained by the above method has data adaptability, which can both filter high noise points and retain micro-disturbance signals, effectively supporting the subsequent dynamic identification of micro-hazard areas.

[0109] The final weak disturbance feature point set D1 is organized into a structured data table, including the following fields: sensor number or monitoring point number; spatial coordinates (x, y, z); periodic index P value; RMS disturbance amplitude; spatial grid number; and time series integrity score.

[0110] A response evolution curve is a continuous function describing the response trend of a monitoring point by mathematically fitting data over an operating period. This process specifically includes the following steps:

[0111] For each data point in D1, the monitoring value within its complete operating cycle is extracted to form the original time series S(t). To enhance data stability, S(t) is first filled with cubic spline interpolation, with the interpolation accuracy controlled within half of the original sampling frequency to ensure smooth and continuous data.

[0112] This invention preferably uses a multinomial regression model to fit S(t), and the fitting function is expressed as follows: ;in, to denoted by , where is the fitting coefficient and is the fitting order. The order 'n' is selected using the Bayesian Information Criterion (BIC), which automatically optimizes the fit error and model complexity to avoid overfitting and underfitting.

[0113] To take into account both global trends and local perturbation changes, a sliding window local fitting mechanism is introduced on the basis of the initial fitting. That is, the curve slides on the original curve with a fixed time window length (such as 24 hours) and performs local linear regression for each window segment to enhance the sensitivity to changes in the response slope.

[0114] The final fitted response curve S*(t) needs to be normalized to linearly map its numerical range to the [0,1] interval. The standardized curve is more conducive to comparison between different points and eliminates sensor calibration differences.

[0115] After obtaining S*(t), mathematical features are extracted to construct a set of hidden danger feature parameters. Specifically, the extracted features include the following four categories:

[0116] Calculate the first derivative of S*(t), representing the rate of change per unit time. Extract its maximum rate of change k_max and the time t1 at which it occurs. This parameter reflects the degree of abrupt change in the response; the steeper the change, the more significant the risk response.

[0117] Periodic analysis of S*(t) using the autocorrelation function method identifies the time lag τ corresponding to the most significant autocorrelation peak, denoted as the main fluctuation period T_p. This parameter reflects the periodic behavior in the response and can be used to identify periodic disturbances driven by external boundary conditions such as air temperature and water level.

[0118] The time interval between the first significant change point t0 of the curve and its maximum response point t_max is defined as the delay time Δt = t_max − t0. The delay time can be used to assess the degree of hysteresis in the response of an engineering structure to external changes. If Δt is long, it may indicate an anomaly in stress transmission.

[0119] Statistical analysis is performed on all maxima and minima in the normalized curve, and indicators such as the number of extreme points N_e, the average extreme value spacing Δx_avg, and the extreme value amplitude difference A_diff_avg are recorded to reflect the degree of instability of the structural response.

[0120] Each data point will correspond to a set of five-dimensional feature vectors, containing a total of five scalar values ​​from the four types of features mentioned above, which are used to construct its potential behavior pattern.

[0121] After summing the feature vectors of all data points in D1, the original set of potential hazard parameters P_raw is formed. To ensure the consistency of subsequent model inputs, this set needs to be standardized and dimensionality reduced. The specific steps are as follows:

[0122] For each feature dimension in P_raw, the Z-score method is used for standardization. The standardization formula is: Z = (X− μ) ÷ σ; where X is the original feature value, μ is the mean of that dimension, and σ is the standard deviation. After standardization, all features have a standard distribution with zero mean and unit variance.

[0123] The processed parameter set is uniformly named the hazard feature parameter set P, which is used for coupling mapping with the spatial model in subsequent steps. This parameter set has a matrix structure, with the number of rows equal to the number of data points in D1, and the number of columns being the feature dimensions retained after standardization.

[0124] In the data subset D1, a set of multi-dimensional feature vectors representing the response behavior is extracted from each data point, forming a hazard feature parameter set P. To make this parameter set locatable, the parameter vector of each data point needs to be associated with its physical acquisition location.

[0125] When collecting monitoring data, each sensor or image acquisition point has clearly defined spatial coordinates (x, y, z), which originate from the engineering layout drawings or BIM model. When forming the feature parameter set P, it needs to be mapped one-to-one with the corresponding coordinates to construct a spatially labeled feature point set Ps, where each element includes:

[0126] Characteristic parameters include point number, spatial location (x, y, z), slope of change, fluctuation period, response delay time, number of local extrema, and difference in extremum amplitude.

[0127] Since digital twin models are typically built on engineering coordinate systems (such as geographic coordinates or structural design coordinates), the collected coordinates may have offsets or differences in accuracy, thus requiring coordinate unification. It is preferable to use a rigid transformation matrix or the ICP (Iterative Closest Point) algorithm to align the original coordinate system to the twin model's coordinate system.

[0128] After coordinate mapping is completed, each feature point is written to the database as an attribute supplement to the twin model mesh (or structural unit). If the twin model uses a structural mesh as its basic unit, point data can be assigned to the corresponding unit using a spatial indexing algorithm (such as a KD-tree). Ultimately, each mesh unit in the model can contain the following information:

[0129] Grid ID, number of feature points, mean of principal feature values, maximum slope, minimum period, risk weight value, etc.

[0130] Based on the main feature parameters (such as the slope of change or the delay time), point data or grid data are rendered into a 3D model using color gradients to form a preliminary spatial heat map. The higher the risk index value, the darker the color of the area, providing an intuitive basis for subsequent regional clustering analysis.

[0131] After feature points are spatially injected, further clustering analysis is needed to examine their distribution and feature consistency in order to identify local regions with trend aggregation. This process mainly includes the following steps.

[0132] For each feature point mapped into the model, a fixed-radius spherical neighborhood method is used, defining its neighborhood radius R (e.g., 5 meters). Other feature points are searched within this sphere to form a neighborhood set. This operation can be accelerated using KD-trees or octrees.

[0133] For each feature point and its neighborhood data, calculate the following two metrics:

[0134] High-risk point density ρ in the neighborhood: defined as the proportion of points that satisfy a certain feature (such as a slope of change greater than the threshold k0);

[0135] Intra-neighborhood feature consistency σ: defined as the standard deviation of a certain principal feature (such as fluctuation period) of all points in the neighborhood.

[0136] If the density ρ in a certain neighborhood is higher than the threshold ρ0 (e.g., 0.7) and the standard deviation σ is lower than σ0 (e.g., 0.1), then the neighborhood is considered a candidate for anomaly region with "consistent characteristic trends".

[0137] Input the set of all points that meet the above conditions into a clustering algorithm for spatial aggregation. The density-based clustering method DBSCAN is preferred, with a minimum number of points minPts (e.g., 5) and a maximum radius ε (e.g., 3 meters) set to extract dense, highly consistent regions in the space.

[0138] Each clustering result represents a potential cluster of hidden danger trends, and is assigned a unique identifier ID. Its boundary point coordinates, number of feature points, and statistical mean of main features are extracted.

[0139] The identified cluster areas need to be clearly marked in the digital twin model and used as input for subsequent simulation calculations and inspection command scheduling.

[0140] Each cluster region is assigned a unique number (such as R_001, R_002, etc.) and accompanied by the following information: region voxel boundary coordinates, number of feature points contained, maximum risk index value, average delay time, etc.

[0141] A new "potential anomaly region identifier" field is added to the structural unit data structure of the twin model, and the region number is bound to its covering unit to form a spatial-attribute-risk multidimensional mapping structure.

[0142] Store complete attribute information for each abnormal region in the model database or nested structure.

[0143] By using 3D visualization tools (such as Unity or Cesium), the marked areas can be presented in the twin model as semi-transparent red areas or dynamic thermal color blocks, enabling engineers to intuitively identify and interact with risk areas.

[0144] Within each identified potential anomaly region, multiple feature points exist, all of which possess historical monitoring time-series data. To reflect the evolutionary trend of the overall structural response in this region, curve modeling of these data is necessary.

[0145] First, extract time-series data of all feature points within the potential anomaly region from the database, including at least:

[0146] Structural strain time series: reflects the deformation of a structure under external loads or internal forces;

[0147] Pore ​​water pressure time series: reflects the state of seepage or pressure changes inside the dam body;

[0148] Microseismic response time series: contains vibration amplitude and frequency information, reflecting subtle structural disturbances;

[0149] Temperature change time series: used to eliminate erroneous responses caused by thermal expansion and contraction;

[0150] Inspection image analysis index sequence: such as the change of crack length and width over time.

[0151] To improve the accuracy of trend identification, the original time series needs to be processed as follows:

[0152] Use polynomial trend lines or moving averages to remove long-term linear drift.

[0153] Synchronize data with different sampling frequencies and unify the sampling interval (e.g., once per hour).

[0154] Missing data were filled using cubic spline interpolation, and outliers were replaced using the median.

[0155] Normalize physical quantities to ensure dimensional uniformity and facilitate characteristic comparison.

[0156] The processed data is smoothed using local weighted regression (LOWESS) to obtain a continuous response curve. Each physical quantity forms a structural change trend curve, which together constitute the trend curve set for that region.

[0157] Construct correlation matrices between trend curves, such as calculating the Pearson correlation coefficient between strain curves and pore pressure curves, and the mutual information index between microseismic frequency and crack propagation, to identify which physical changes are driven by coupling and to determine whether the damage is induced by a single factor or caused by the superposition of multiple factors.

[0158] The evolution of the trend curve is determined not only by the state of the structure itself, but also by the influence of external environmental boundary conditions. Therefore, this invention further incorporates environmental data to invert the micro-damage type and distribution direction.

[0159] Boundary conditions include, but are not limited to: daily average temperature and temperature gradient; reservoir water level change rate; daily rainfall and cumulative rainfall; groundwater level and reservoir water pressure; and seasonal geothermal variation trends. All boundary data are aligned to the time axis of the trend curve at a uniform time step.

[0160] The response coupling relationship between the trend curve and boundary conditions is analyzed using the following indicators:

[0161] Pearson correlation coefficient (r): used to determine whether the relationship is linear;

[0162] Mutual information value (MI): used to determine the strength of nonlinear driving relationships;

[0163] Time delay analysis value: used to identify hysteresis response, such as strain delay caused by temperature change.

[0164] If the structural response curve is significantly correlated with a specific boundary condition (e.g., r is greater than 0.6 or MI value is greater than 0.3), it can be inferred that the response in that region is closely related to that boundary factor.

[0165] Based on the existing micro-damage database, load typical micro-damage response templates, including but not limited to:

[0166] Opening crack mode: characterized by rapid increase in strain and decrease in pore pressure;

[0167] Shear slip mode: characterized by increased microseismic frequency and synchronous strain changes in multiple directions;

[0168] The seepage scour mode is characterized by frequent pore pressure fluctuations and enhanced micro-vibration amplitude.

[0169] Structural relaxation mode: characterized by delayed response and increased low-frequency amplitude.

[0170] By matching the trend curve set with typical templates using the minimum mean square error or maximum correlation coefficient method, the micro-damage type that best fits the region can be obtained.

[0171] Based on the spatial distribution direction of the trend response (such as along the dam axis, perpendicular to the dam surface, and slope direction) combined with the direction of the structural strain tensor, the direction of defect action and evolution is determined, forming spatial descriptive attributes.

[0172] The estimated micro-damage type, its spatial location, and direction of action are structured and organized to generate the inverted defect feature set Q.

[0173] Each defective unit contains the following fields:

[0174] Defect types (opening, shearing, scouring, etc.);

[0175] Spatial location (x, y, z);

[0176] The main direction of the defect (such as the negative X-axis direction, the 45° slope direction, etc.);

[0177] Damage intensity indicators (such as maximum strain increment, microseismic energy, etc.);

[0178] Description of damage trend (rising, stable, drastic fluctuations, etc.).

[0179] All defect information is stored in the form of a structured list or a JSON object.

[0180] All defect points are combined into an inverted defect feature set Q, which serves as an important input for the initial defect boundary conditions in the subsequent local simulation model, and also as incremental information for the dynamic evolution and update of the twin model.

[0181] The inverted defect feature set Q is obtained from the previous steps and records the location coordinates, type, direction of action, risk intensity, and trend parameters of each potential micro-damage.

[0182] For each defect record, its spatial location (x, y, z) is used as a geometric embedding point. Combined with its damage type, a corresponding structural perturbation entity is constructed. An example is shown below:

[0183] Crack-type defects are modeled as "zero-thickness interfaces" by defining contact surface elements;

[0184] For material softening defect regions, parameter perturbation can be achieved by reducing the elastic modulus or introducing damage constitutive relations.

[0185] Pore ​​pressure anomaly defects are addressed by setting up local high-pressure sources or permeability abrupt change zones in the seepage field.

[0186] A local simulation sub-region is constructed around each defect point, typically with a spatial extent 3 to 5 times the radius of the defect's influence. Tetrahedral or hexahedral finite element methods are preferred for mesh generation, with anisotropic mesh refinement based on the principal direction of the defect to enhance the simulation accuracy of stress concentration zones.

[0187] The currently measured temperature field, water pressure field, and structural load are used as the initial field inputs, while boundary conditions are set for a future preset period (e.g., 7 days, 30 days, or 90 days), including:

[0188] External water level change curve;

[0189] Ambient temperature prediction curve;

[0190] Simulation of load changes (such as upstream water storage, flood discharge, seasonal deformation, etc.).

[0191] Boundary condition data can be provided through historical model predictions, engineering operation plans, or weather forecast models.

[0192] Based on the completed local model, a complete set of multi-field coupled control equations is established and solved using time-stepping.

[0193] This invention employs the following three sets of governing equations:

[0194] The heat conduction equation: The governing equation for the change of temperature T with time t is: specific heat multiplied by density multiplied by the first derivative of temperature with respect to time equals thermal conductivity multiplied by the second derivative of temperature with respect to spatial coordinates.

[0195] The seepage control equation is as follows: the change of pore pressure P follows Darcy's law and the continuity equation, which can be expressed as: the change of pore pressure with time is equal to the product of the permeability coefficient and the pressure gradient, and is affected by structural deformation.

[0196] Stress equilibrium equation: The structural stress σ follows the force equilibrium condition, that is: the derivative of the internal stress tensor with respect to spatial coordinates is equal to the body force density. Combined with constitutive relations (such as the Drucker-Prager model), it reflects the nonlinear response of the material.

[0197] The three fields mentioned above interact through coupling terms. For example, temperature affects material expansion, thus inducing stress changes; the osmotic pressure field affects the effective stress of the structure; and stress deformation, in turn, affects porosity, thereby affecting the seepage path. The coupling strategy can be implemented using "stepwise iterative coupling" or "strongly coupled element method". This invention preferably adopts the element strong coupling method to reduce error propagation.

[0198] Spatial discretization is performed using the finite element method, with a backward Euler scheme used for the time dimension. The time step Δt (e.g., 1 hour) is set, and the total simulation time T (e.g., 720 hours) requires 720 iterations. After each solution step, the principal values ​​of the stress tensor for each structural element within the region are recorded. And form time series data.

[0199] After the simulation is completed, the stress-time curves at key locations in the potential anomaly region are extracted, i.e., the stress response evolution curves.

[0200] The original output may contain discrete oscillations and fluctuations, which need to be processed using the following method:

[0201] The sliding window averaging method was used (window width was 5 to 10 time steps).

[0202] Perform first-order low-pass filtering to preserve the main trend;

[0203] Finally, a continuous and highly readable stress evolution curve S(t) is obtained.

[0204] Extract the following indicators from S(t):

[0205] Maximum principal stress σ_max;

[0206] Average growth rate (stress increase per unit time);

[0207] Stress growth acceleration (i.e., the slope of the slope);

[0208] Stress trend fitting types (linear, power function, exponential growth, etc.);

[0209] The trend type can be selected based on the fitted residuals. Commonly used fitting functions are: Where a, b, and c are fitting parameters, and b reflects the degree of acceleration.

[0210] The curves are embedded into the digital twin system interface in the form of a two-dimensional chart, while areas of accelerated trend are marked on the three-dimensional model with risk color labels.

[0211] All output data is written to the database to support subsequent risk score calculation and dynamic early warning analysis.

[0212] After completing the local multi-field simulation, the stress principal value time series curve S(t) of key points in the potential anomaly region is obtained, which is expressed as: ,in Let be the maximum principal stress, t be time, and n be the total number of time steps. To comprehensively evaluate the changing trend and stability of the curve, this invention extracts key features from the following three dimensions:

[0213] The first derivative of the curve S(t), i.e., the rate of stress change per unit time, is calculated using the formula: R = maximum stress growth rate ÷ average stress change rate; where the maximum growth rate is the maximum value among all Δσ / Δt, and the average change rate is the total stress increment divided by the total time. If the R value is significantly greater than 1, it indicates that the curve shows an accelerating upward trend, and the structural risk is amplified.

[0214] The mean difference between adjacent windows is calculated using the sliding window method (window length W, usually set to 5 to 10 time steps). If a certain difference Δσ is greater than a set threshold value (such as twice the baseline fluctuation amplitude), it is identified as a "step event".

[0215] After accumulating all step events, the step intensity index S is defined as: S = total number of step events × average step amplitude; this index reflects whether there are sudden changes in the structural response, which may be caused by boundary disturbances, structural failure or critical instability.

[0216] Perform a Fast Fourier Transform (FFT) or a continuous wavelet transform on S(t) to obtain its frequency component spectrum. After identifying frequency abrupt changes (such as dominant frequency jumps or bandwidth expansion), define the frequency perturbation factor as: F_var = energy percentage of the frequency abrupt change segment × dominant frequency offset amplitude; this factor is used to reflect whether abnormal frequency domain perturbations have occurred in the structural response, and is usually associated with physical mechanisms such as material softening and boundary degradation.

[0217] After obtaining the three basic indicators R, S, and F_var, this invention uses a weighted linear combination to construct the risk score value V: ;in: Let be the weighting coefficient, satisfying The weights can be set based on engineering experience or optimized through training with historical cases. An example is set as follows: V represents the final risk score, which is determined based on the normalization method.

[0218] By retrospectively analyzing typical engineering accident cases and daily operational data, the corresponding V-value distribution was statistically analyzed, and the scoring space was automatically divided using K-means clustering or piecewise regression. The following thresholds were ultimately set:

[0219] T1: The dividing line between low and medium risk, a value of 0.25 is recommended;

[0220] T2: The dividing line between medium and high risk, it is recommended to take 0.50;

[0221] T3: The dividing line between high risk and extremely high risk; a value of 0.75 is recommended.

[0222] Users can also adjust the model based on their actual project experience to achieve adaptability in different projects.

[0223] When the calculated risk score V meets the following conditions, the corresponding risk level L will be output:

[0224] If V ≤ T1, then it is considered "low risk";

[0225] If T1 < V ≤ T2, then it is considered "medium risk";

[0226] If T2 < V ≤ T3, then it is "high risk";

[0227] If V > T3, then it is "extremely high risk".

[0228] Each risk level will be accompanied by trend information as a label, including:

[0229] An upward trend (R continues to increase);

[0230] Fluctuation trend (S appears frequently);

[0231] Frequency perturbation increased (F_var increased abruptly);

[0232] Multiple factors are superimposed (R, S, and F_var are all high).

[0233] Risk levels are mapped onto the surface of the digital twin 3D model using a color gradient. For example:

[0234] Green indicates low risk;

[0235] Yellow indicates medium risk;

[0236] Orange indicates high risk;

[0237] Red or purple indicates extremely high risk.

[0238] Color mapping follows a linear or exponential proportional relationship, making it easy to intuitively judge the structural state.

[0239] Trend labels (such as "rapid growth", "frequency disorder", "multi-factor overlay") are displayed on the surface of the corresponding area or in a floating window using arrows, labels or dynamic icons.

[0240] Each risk score V, level L, and trend description is written into the structural unit attribute table of the twin model for model state version management and evolution path comparison, and can trigger the following operations:

[0241] Risk alerts are pushed out; inspection or testing suggestions are automatically generated; emergency response or simulation prediction tasks are initiated.

[0242] Example 2, please refer to Figure 2 As shown in this embodiment, the water conservancy project monitoring and analysis system based on digital twins includes:

[0243] Data acquisition module: Collects multi-source monitoring data during the operation cycle of the target water conservancy project, and synchronizes and aligns the data according to spatial location and time order to construct a standardized monitoring dataset D;

[0244] Data subset construction module: Based on the historical change patterns in dataset D, identify data points with periodic fluctuation characteristics, and extract a set of data subsets D1 with weak perturbation characteristics by combining their spatial distribution density;

[0245] Hazard feature extraction module: Perform multi-parameter joint fitting on data subset D1, establish the response evolution curve of each data point, and extract the curve's slope, fluctuation period, delay time, and local extreme value feature parameters to form a hazard feature parameter set;

[0246] Potential anomaly region marking module: Maps the set of hazard characteristic parameters to the spatial coordinate system of the digital twin model, filters out local areas with abnormal evolution trends, and marks them as potential anomaly regions in the digital twin model;

[0247] Inversion Defect Module: Based on the historical time series data corresponding to the potential anomaly region, construct the local structural change trend curve of the potential anomaly region, and inversely deduce the possible micro-damage types and their distribution directions by combining the current environmental boundary conditions, thus obtaining the inversion defect feature set;

[0248] Coupled calculation module: Input the inverted defect feature set into the digital twin model, perform local thermal-permeability-mechanical coupled calculation, and obtain the stress response evolution curve of the potential anomaly area in the future preset period;

[0249] Risk level classification module: Based on the growth rate, step characteristics and frequency mutation in the stress response curve, calculate the risk score value V of the potential anomaly area. If V exceeds the set threshold T, output the risk level of the potential anomaly area and the corresponding evolution trend information.

[0250] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A water conservancy project supervision analysis method based on digital twinning, characterized in that: The method comprises the following steps: Collecting multi-source monitoring data of the target water conservancy project in the operation period, and synchronously aligning the data according to spatial position and time sequence to construct a standardized monitoring data set D; Based on the historical change law in the data set D, data points with periodic fluctuation characteristics are identified, and a data subset D1 with weak disturbance characteristics is extracted by combining the spatial distribution density; The data subset D1 is fitted by multiple parameters to establish the response evolution curve of each data point, and the change slope, fluctuation period, delay time and local extreme value characteristic parameters of the curve are extracted to form a hidden danger characteristic parameter set; The hidden danger characteristic parameter set is mapped to the spatial coordinate system of the digital twin model, the local area with abnormal evolution trend is screened out, and the potential abnormal area is marked in the digital twin model; Based on the corresponding historical time series data in the potential abnormal area, the local structure change trend curve of the potential abnormal area is constructed, and the possible micro-damage type and its distribution direction are inversely deduced combined with the current environmental boundary conditions to obtain the inversion defect characteristic set; The construction of the local structure change trend curve of the potential abnormal area based on the corresponding historical time series data in the potential abnormal area comprises the following steps: extracting the time series data of all feature points in the potential abnormal area in the historical operation period, including strain, pore water pressure, microseismic response amplitude and frequency parameters; after the long-term stable drift is eliminated by detrending processing of the time series of each type of physical quantity, the sliding average and local weighted regression method is used to construct the smooth response curve; the time change gradient, change rate and key turning point of each type of physical quantity are extracted by curve fitting to construct a set of multivariate structure response trend curves; the multiple physical field response curves are aligned according to the time axis, and an association matrix is established to identify the cooperative relationship between the changes of different physical quantities; The correlation analysis between the local structure change trend curve and the current boundary conditions is carried out, the influence intensity is quantified by using the Pearson correlation coefficient and the mutual information index; combined with the multiple typical response modes in the preset micro-damage mechanism library, the current trend curve is matched, including: according to the existing micro-damage database, loading the typical micro-damage response template, matching the trend curve group with the typical micro-damage response template by the least mean square error or the maximum correlation coefficient method, obtaining the micro-damage type that meets the area, and recording the occurrence direction and spatial offset position; The inversion defect characteristic set is input into the digital twin model, the local heat conduction-seepage-stress coupling calculation is performed, and the stress response evolution curve of the potential abnormal area in the future preset period is obtained; According to the growth rate, step feature and frequency mutation in the stress response curve, the risk score value V of the potential abnormal area is calculated, and if V exceeds the set threshold T, the risk level and the corresponding evolution trend information of the potential abnormal area are output.

2. The hydraulic engineering supervision analysis method based on digital twinning according to claim 1, characterized in that: The multi-source monitoring data includes dam body surface temperature, structural strain, pore water pressure, microseismic amplitude and frequency, and crack features in artificial inspection images.

3. The hydraulic engineering supervision analysis method based on digital twinning according to claim 2, characterized in that: The data points with periodic fluctuation characteristics are identified, including: A time series curve is constructed for each monitoring point in the standardized monitoring data set, and a first derivative thereof is calculated using a sliding time window to extract a local change rate; A frequency domain analysis is performed on the time series curve using a fast Fourier transform to identify a main frequency component and calculate a periodicity index P, the periodicity index P being a ratio of a main frequency amplitude to an average energy; Monitoring points with a periodicity index P greater than a set threshold P0 are determined to be data points with periodic fluctuation characteristics.

4. The hydraulic engineering supervision and analysis method based on digital twinning according to claim 3, characterized in that: The multi-parameter joint fitting of the data subset D1 includes: For each data point in D1, a complete time series thereof is extracted, and missing data is completed using a cubic spline interpolation method; A polynomial regression model is used to preliminarily fit the time series, and the order of the fitting function is determined according to the principle of minimum Bayesian information criterion; On the basis of the regression fitting, a sliding window method is further used to locally optimize the curve shape, extract a local change trend, and construct a complete response evolution curve; The fitted curve is normalized to standardize the response change range of all data points to the interval [0, 1].

5. The hydraulic engineering supervision and analysis method based on digital twinning according to claim 4, characterized in that: The change slope, fluctuation period, delay time and local extreme value characteristic parameters of the curve are extracted, including: The change rate function of the response curve is obtained by calculating the first derivative of the response curve in the whole period, so as to extract the maximum change slope value and the corresponding time position, which are used to reflect the rate characteristics of the response intensity rising stage; Periodic analysis is performed on the normalized curve, and the autocorrelation function method is used to identify the main period, and the period length is recorded as the fluctuation period feature; The time difference between the initial disturbance point and the maximum response point is calculated by a time delay analysis method to extract the delay time parameter; The positions and amplitudes of all local maximum and minimum points in the curve are identified, and the number of extreme points, average spacing and amplitude difference are counted to form a local extreme value feature subset.

6. The hydraulic engineering supervision and analysis method based on digital twinning according to claim 5, characterized in that: The local area with an abnormal evolution trend is screened out, including: The spatial position coordinate information (x, y, z) of each hidden danger feature vector corresponding to the data acquisition stage is bound to construct a feature point set with a spatial label; According to the feature point set, a fixed radius sphere neighborhood method is used to delimit the local neighborhood of each point, and the mean and variance of the risk features in the neighborhood are counted; A risk density function is constructed, and if the high-risk point density in a neighborhood exceeds a set threshold p0 and the standard deviation of the feature values in the region is lower than a threshold s0, it is judged to be an abnormal region with consistent trend; A clustering algorithm is used to spatially cluster the point set that meets the conditions to extract multiple local risk regions with consistent evolution; Based on the clustering results, the boundary voxel set of each region is extracted to form a clear spatial boundary of the abnormal region.

7. The hydraulic engineering supervision and analysis method based on digital twinning according to claim 1, characterized in that: The local heat conduction-seepage-stress coupling calculation is performed, including: The spatial boundary range and time scale of the coupling solution region are set, and the boundary range includes the defect point and the structural units within a range of not less than 3 times the influence radius around the defect point. The temperature, water pressure, and structural load at the current time are introduced as initial boundary conditions, and the future meteorological and operation prediction data for a preset period are combined to construct the time step driven boundary; The numerical solution of the heat conduction-seepage-stress coupling equation set is realized based on the finite element method or the finite volume method; At each time step, the principal value change of the stress tensor of the local defect area is extracted to form a stress response evolution curve, which is used to represent the structural evolution trend.

8. A water conservancy project supervision analysis system based on digital twinning, for implementing the water conservancy project supervision analysis method based on digital twinning according to any one of claims 1-7, characterized in that: Comprise: Data acquisition module: collect multi-source monitoring data of the target water conservancy project during the operation period, and synchronize and align them according to spatial position and time sequence to construct a standardized monitoring data set D; Data subset establishment module: based on the historical change law in the data set D, identify data points with periodic fluctuation characteristics, and extract a data subset D1 with weak disturbance characteristics based on the spatial distribution density; Hidden danger feature extraction module: multi-parameter joint fitting is performed on the data subset D1 to establish the response evolution curve of each data point, and the change slope, fluctuation period, delay time, and local extreme feature parameters of the curve are extracted to form a hidden danger feature parameter set; Potential abnormal area marking module: map the hidden danger feature parameter set to the spatial coordinate system of the digital twin model, filter out local areas with abnormal evolution trend, and mark them as potential abnormal areas in the digital twin model; Inversion defect module: based on the corresponding historical time series data in the potential abnormal area, construct the local structural change trend curve of the potential abnormal area, and inversely deduce the possible types and distribution directions of micro-damage based on the current environmental boundary conditions to obtain an inversion defect feature set; Coupling calculation module: input the inversion defect feature set into the digital twin model to perform local thermal-seepage-force coupling calculation and obtain the stress response evolution curve of the potential abnormal area in the future preset period; Risk level division module: according to the growth rate, step feature, and frequency mutation in the stress response curve, calculate the risk score value V of the potential abnormal area, if V exceeds the set threshold T, output the risk level and corresponding evolution trend information of the potential abnormal area.

Images