Hydraulic engineering supervision and analysis system and method based on digital twinning

Abstract

The invention discloses a hydraulic engineering supervision and analysis system and method based on digital twinning, and belongs to the technical field of hydraulic engineering operation monitoring and management, and the method comprises the steps: recognizing data points with periodic fluctuation characteristics and sparse spatial distribution, extracting a weak disturbance data subset, and building a response evolution curve; identifying a potential abnormal region through neighborhood statistics and clustering; for the abnormal region, inverting possibly existing micro-damage types and distribution directions based on historical time sequences and environmental boundary conditions to obtain an inversion defect feature set, and inputting the inversion defect feature set into a twinborn model to execute local heat-seepage-force coupling calculation to obtain a stress response evolution curve in a future preset period; calculating a risk score value according to the growth rate, step characteristics and frequency abrupt change of the curve, and outputting risk level and evolution trend information when the score exceeds a set threshold value; according to the invention, early recognition and dynamic evolution analysis of tiny hidden dangers can be realized, and the intelligent and refined level of hydraulic engineering safety supervision is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of water conservancy project operation monitoring and management, in particular to a water conservancy project supervision analysis system and method based on digital twinning. BACKGROUND

[0002] Water conservancy projects such as dams, sluices, reservoirs, etc., as important infrastructure for controlling water flow, water storage, flood control and disaster reduction, and irrigation water supply, their safe operation has a significant impact on regional economic development and the safety of life and property of residents. Traditional water conservancy project supervision mainly relies on on-site inspection, manual inspection or monitoring based on limited sensor data, which has problems such as poor timeliness, low perception dimension, and lagging fault response, etc., and it is difficult to meet the requirements of fine control and early warning of equipment operating conditions under complex working conditions.

[0003] Especially in special working conditions such as seasonal alternation and hydrological mutation, water conservancy projects may have local abnormal hidden dangers such as piping, leakage, and dam body cracks. Since these hidden dangers usually have very small parameter changes at the initial stage, which are often below the equipment sensing threshold, and have irregularity and instantaneity in spatial distribution, traditional methods cannot achieve early identification. For example, some soil dams have the potential risk of crack expansion in the dam shoulder section after the dry season to the early flood season, and such small hidden dangers are difficult to form effective features in sensor response, becoming a "blind area" that needs to be solved in water conservancy project safety supervision.

[0004] Digital twinning technology, as an advanced simulation method that integrates sensor data, historical working conditions, engineering modeling and AI analysis, provides full life cycle modeling, prediction and feedback capabilities for complex engineering structures. However, the current digital twinning system in water conservancy projects mainly stays at the level of three-dimensional modeling and state playback, lacking identification and response mechanisms for abnormal features, and cannot accurately model and predict small but critical safety hazards, resulting in a "perception gap" between monitoring data and actual safety status. SUMMARY

[0005] The purpose of the present application is to provide a water conservancy project supervision analysis system and method based on digital twinning to solve the problems in the background art.

[0006] In order to achieve the above purpose, the present application provides the following technical solution: a water conservancy project supervision analysis method based on digital twinning, comprising: Collecting multi-source monitoring data of the target water conservancy project in the operation cycle, and synchronously aligning them 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, identifying data points with periodic fluctuation characteristics, and combining their spatial distribution density to extract a data subset D1 with weak disturbance characteristics; The multi-parameter joint fitting is performed on the data subset D1, an evolution curve of response of each data point is established, and a change slope, a fluctuation period, a delay time and a local extreme characteristic parameter of the curve are extracted to form a hidden danger characteristic parameter set; The hidden danger characteristic parameter set is mapped into a spatial coordinate system of the digital twin model, a local area with an abnormal evolution trend is screened out, and is marked as a potential abnormal area in the digital twin model; Based on the corresponding historical time series data in the potential abnormal area, a local structural change trend curve of the potential abnormal area is constructed, and a possible micro-damage type and its distribution direction are backstepped combined with a current environmental boundary condition to obtain an inversion defect characteristic set; The inversion defect characteristic set is input into the digital twin model, a local thermal-infiltration-force coupling calculation is performed, and a stress response evolution curve of the potential abnormal area in a future preset period is obtained; According to a growth rate, a step feature and a frequency mutation in the stress response curve, a risk score value V of the potential abnormal area is calculated, and if V exceeds a set threshold T, a risk level and corresponding evolution trend information of the potential abnormal area are output.

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

[0008] Preferably, 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 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, a main frequency component is identified, and a periodicity index P is calculated, the periodicity index P being a ratio of a main frequency amplitude to an average energy; The monitoring points with a periodicity index P greater than a set threshold P0 are determined as data points with periodic fluctuation characteristics.

[0009] Preferably, the multi-parameter joint fitting on 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 a degree of the fitting function is determined according to a minimum principle of a Bayesian information criterion; On the basis of the regression fitting, a sliding window method is used to further locally optimize a curve shape, a local change trend is extracted, and a complete response evolution curve is constructed; The fitted curve is normalized to standardize a response change range of all data points to an interval [0, 1].

[0010] Preferably, the change slope, fluctuation period, delay time and local extreme characteristic parameters of the extraction curve include: By calculating the first derivative of the response curve in the whole cycle, the change rate function is obtained, and the maximum change slope value and its corresponding time position are extracted, which are used to reflect the rate characteristics of the rising stage of response intensity; Periodic analysis is performed on the normalized curve, the autocorrelation function method is used to identify the main period, and the period length is recorded as the fluctuation period characteristic; The time difference between the initial disturbance point and the maximum response point is calculated by the time delay analysis method, and the delay time parameter is extracted; Identify the position and amplitude of all local maximum and minimum points in the curve, and count the number of extreme points, average spacing and amplitude difference to form a local extreme characteristic subset.

[0011] Preferably, the local area with abnormal evolution trend is screened out, including: Bind the spatial position coordinate information (x, y, z) of each hidden danger characteristic vector corresponding to the data acquisition stage to construct a feature point set with spatial tags; According to the feature point set, a local neighborhood of each point is determined by using the fixed radius sphere neighborhood method, and the mean and variance of the risk characteristics in the neighborhood are counted; Construct a risk density function. If the high-risk point density in a neighborhood exceeds the set threshold ρ0, and the standard deviation of the eigenvalues in the region is lower than the threshold σ0, it is judged as an abnormal region with consistent trend; Using clustering algorithm to cluster the point set that meets the condition, extract multiple local risk areas with consistent evolution; Based on the clustering results, extract the boundary voxel set of each region to form a clear spatial boundary of the abnormal region.

[0012] Preferably, the local structural change trend curve of the potential abnormal area is constructed based on the corresponding historical time series data in the potential abnormal area, including: Extract the time series data of all feature points in the potential abnormal area in the historical operation cycle, including strain, pore water pressure, microseismic response amplitude and frequency parameters; After detrending the time series of each type of physical quantity, eliminating long-term steady-state drift, and using moving average and local weighted regression method to construct smooth response curve; Extract the time change gradient, change rate and key turning point of each type of physical quantity by curve fitting to construct a set of multivariate structural response trend curves; Align the response curves of multiple physical fields according to the time axis, and establish a correlation matrix to identify the cooperative relationship between the changes of different physical quantities.

[0013] Preferably, the possible micro-damage type and its distribution direction are deduced by combining the current environmental boundary conditions, including: Obtaining the current environmental boundary conditions corresponding to the abnormal area, including air temperature, water level, reservoir water pressure change rate, rainfall and seasonal ground temperature gradient; Correlation analysis is performed on the local structure change trend curve and the current boundary condition, and the influence intensity is quantified by using Pearson correlation coefficient and mutual information index; Combining the multiple typical response modes in the preset micro-damage mechanism library, the current trend curve is matched with the mode; The most possible micro-damage type is selected by the minimum residual matching principle, and its occurrence direction and spatial offset position are recorded.

[0014] Preferably, the local thermal-seepage-force 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 structure unit with a peripheral range not less than 3 times the influence radius; The temperature, water pressure and structure load at the current time are imported as initial boundary conditions, and the time step driving boundary is constructed by combining the future preset period of meteorological and operation prediction data; 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; The principal value change of the stress tensor of the local defect area at each time step is extracted to form a stress response evolution curve, which is used to represent the structure evolution trend.

[0015] The application also provides a water conservancy project supervision analysis system based on digital twinning, comprising: A data acquisition module: acquiring multi-source monitoring data in the operation period of the target water conservancy project, and synchronously aligning the data according to spatial position and time sequence to construct a standardized monitoring data set D; A data subset establishment module: based on the historical change law in the data set D, identifying data points with periodic fluctuation characteristics, and combining the spatial distribution density to extract a data subset D1 with weak disturbance characteristics; A 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 value feature parameters of the curve are extracted to form a hidden danger feature parameter set; A potential abnormal area marking module: mapping the hidden danger feature parameter set to the spatial coordinate system of the digital twinning model, screening out the local area with abnormal evolution trend, and marking it as a potential abnormal area in the digital twinning model; Inversion defect module: based on the corresponding historical time series data in the potential abnormal area, a local structural change trend curve of the potential abnormal area is constructed, and the possible micro-damage type and its distribution direction are back calculated in combination with the current environmental boundary conditions, to obtain an inversion defect feature set; Coupling calculation module: inputting the inversion defect feature set into the digital twin model, performing local thermal-infiltration-force coupling calculation, and obtaining a stress response evolution curve of the potential abnormal area in a future preset period; Risk level division module: according to the growth rate, step feature and frequency mutation in the stress response curve, calculating the risk score value V of the potential abnormal area, and if V exceeds the set threshold T, outputting the risk level and corresponding evolution trend information of the potential abnormal area.

[0016] In the above technical solution, the technical effects and advantages provided by the present application are as follows: 1. The present application realizes quantitative calculation of the structural risk of the potential abnormal area by multi-dimensionally extracting the growth rate, step feature and frequency disturbance factor of the stress response curve, and constructing a weighted risk score model. Compared with the traditional monitoring method relying on a single index (such as stress absolute value or strain rate), the present application can identify the risk trend in advance when the hidden danger is still in the micro-evolution stage, effectively avoiding the missed judgment problem caused by insignificant response, thereby significantly improving the sensitivity and foresight of water conservancy engineering safety supervision.

[0017] 2. The present application constructs a risk score threshold system based on historical operation cases, and realizes dynamic output of risk level and trend information in combination with the digital twin model, forming a closed-loop mechanism of “calculation-determination-visualization-feedback”. This mechanism not only improves the intuitiveness and operability of risk expression, but also provides data-based decision-making basis for inspection and dispatch, emergency disposal and structure reinforcement, which can greatly improve the intelligent and fine level of water conservancy engineering operation safety management. BRIEF DESCRIPTION OF DRAWINGS

[0018] In order to more clearly illustrate the technical solutions in the embodiments or prior art, the drawings needed in the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments described in the present application, and other drawings can also be obtained by those skilled in the art based on these drawings.

[0019] Figure 1 The method mind map of the present application.

[0020] Figure 2 The system module flow chart of the present application. DETAILED DESCRIPTION

[0021] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are some but not all of the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the protection scope of the present application.

[0022] Embodiment 1, please refer to Figure 1 The water conservancy project supervision analysis method based on digital twinning described in the embodiment includes: Collecting multi-source monitoring data in the operation period of the target water conservancy project, and synchronously aligning the monitoring data according to spatial positions and time sequences to construct a standardized monitoring data set D; Based on historical change rules in the data set D, data points with periodic fluctuation characteristics are identified, and a data subset D1 with weak disturbance characteristics is extracted in combination with spatial distribution densities; The data subset D1 is subjected to multi-parameter joint fitting to establish a response evolution curve of each data point, and 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 into a spatial coordinate system of a digital twinning model, a local area with an abnormal evolution trend is screened out, and is marked as a potential abnormal area in the digital twinning model; Based on corresponding historical time sequence data in the potential abnormal area, a local structure change trend curve of the potential abnormal area is constructed, and a possible micro-damage type and its distribution direction are backstepped in combination with current environmental boundary conditions to obtain an inversion defect characteristic set; The inversion defect characteristic set is input into the digital twinning model, local thermal-seepage-force coupling calculation is performed, and a stress response evolution curve of the potential abnormal area in a future preset period is obtained; According to a growth rate, a step feature and a frequency mutation in the stress response curve, a risk score value V of the potential abnormal area is calculated, and if the V exceeds a set threshold T, a risk level and corresponding evolution trend information of the potential abnormal area are output.

[0023] The target water conservancy project mainly refers to a control type hydraulic structure mainly of a concrete or earth-rock structure, such as a dam structure, a sluice and a reservoir, and an operation period thereof can be divided according to seasons, hydrological years or structure maintenance periods. Generally, the operation period is 30 days to 365 days, and is determined according to engineering management requirements.

[0024] In the operation period, the collected monitoring data includes but is not limited to the following types: Dam surface temperature data: Obtain the surface temperature of each key part of the dam body through thermocouples or infrared thermometers, and the sampling frequency is generally 1 time per hour; Structural strain data: Obtained by strain gauges or optical fiber sensors embedded in the concrete structure, with a general sampling frequency of 1 time per 10 minutes; Pore water pressure data: Use pore pressure meters inside and outside the dam body or around the foundation to collect pore water pressure reflecting seepage state, with the same sampling frequency as strain data; Microseismic response data: Use seismic monitoring devices to record the amplitude and frequency characteristics of the dam under external disturbance (such as water release, blasting, earthquake, etc.), with a sampling frequency of more than 10 times per second; Artificial inspection image data: Images obtained by on-site artificial inspection or unmanned aerial vehicle inspection, extracting crack length, width and direction information on the structure surface through image processing algorithms to form a structure crack feature set.

[0025] Due to the large difference in the source of the above-mentioned various data and the different sampling intervals, unified processing is required to ensure that the data can be used for evolution calculation under the same model. Therefore, the present application adopts the following steps to standardize the original data: Data preprocessing: including time format unification, outlier rejection and interpolation filling. All sampling time stamps are converted to Coordinated Universal Time (UTC), and missing points are filled by local linear interpolation; outliers obviously beyond the physical range are replaced by median value processing.

[0026] Time series synchronization alignment: Considering that the sampling frequencies of different monitoring quantities are inconsistent, a unified time step Δt is set, for example, 1 hour. Resample all data according to this time step: For data with sampling frequency higher than Δt (such as microseismic), use moving average; For data with sampling frequency lower than Δt (such as artificial image), keep the latest value or perform linear prediction within the corresponding time window.

[0027] Spatial position normalization: Map the spatial positions corresponding to all monitoring data to a unified three-dimensional coordinate system, taking the dam design drawing or BIM model as the reference. The positions of all sensor installation points are represented in the form of "x, y, z", and the crack coordinates of image data are projected onto the structure model surface through image matching.

[0028] The data obtained after processing is organized into a data set D with each data record including: spatial position coordinates (x, y, z); a time stamp t; a temperature T(x, y, z, t); a strain ε(x, y, z, t); a pore pressure P(x, y, z, t); a microseismic intensity A(x, y, z, t), a frequency f(x, y, z, t); and a crack feature vector C(x, y, z, t) containing crack length, width and strike dimensions.

[0029] The data set D has a fixed sampling interval in time, and all data indicators have uniform physical units and dimensions.

[0030] In actual implementation, the data set D is saved in a database storage structure or a multi-dimensional array form, supports indexing and batch reading, and is used for subsequent hazard identification, local analysis and model evolution simulation.

[0031] For various types of continuous time series data (such as temperature, pore pressure, strain, etc.) in the monitoring data set, the present application first performs periodic analysis to identify data points with stable periodic fluctuation behavior. The following steps are used to achieve this stage: For each sensor or sampling point in the monitoring data set D, the complete time series of the sensor or sampling point in the running period is extracted and set as S(t), where t represents time. The local change rate sequence ΔS(t) is obtained by calculating the difference between adjacent time points, i.e. first-order difference, which is used to eliminate the influence of trend items and highlight short periodic behavior.

[0032] The ΔS(t) sequence is subjected to fast Fourier transform (FFT) to obtain its frequency spectrum information. Let the amplitude of the main frequency component be A_max and the total energy of the frequency spectrum be E_total. The periodicity index P of the point is defined as A_max divided by E_total, i.e.: Periodicity index P = main frequency amplitude A_max ÷ total frequency energy E_total.

[0033] E_total is calculated by summing the energy squares of all frequency components of the frequency spectrum.

[0034] A periodicity identification threshold P0 is set, and data points with a periodicity index P greater than P0 are selected as "periodic fluctuation feature points". The value of P0 is generally between 0.25 and 0.45, and is determined according to the statistical distribution of the frequency spectrum of historical data.

[0035] For the data points selected by the periodicity index, it is further checked whether there is a missing segment exceeding the threshold in the time series. Only data points with a time integrity higher than 95% are retained for subsequent analysis. Finally, a set of "periodic fluctuation feature points" is output for subsequent spatial density analysis.

[0036] Periodic feature points may be concentrated in the structure active area, in order to screen out the micro-disturbance point set with spatial sparse distribution characteristics, further three-dimensional space analysis is needed, the steps are as follows: With the design coordinate system of the dam body as the reference, the whole monitoring area is divided into regular cubic grid, the grid length is L, the recommended value range is 2-5 meters. Each periodic feature point is assigned to the corresponding grid unit according to its spatial coordinates (x, y, z).

[0037] The number of feature points contained in each grid is counted and recorded as The average density N_avg of all grid units is calculated, and the standard deviation σ_N is recorded.

[0038] Set the spatial sparse threshold N_threshold=N_avg−σ_N. Only the points in the grid smaller than N_threshold are retained to form the sparse candidate set.

[0039] For each point in the sparse candidate set, the root mean square value of ΔS(t) is calculated, recorded as RMS_i. Set the disturbance amplitude threshold A0, only the data points with RMS_i less than A0 are retained to form the final weak disturbance data subset D1.

[0040] This process ensures that the data points in D1 have "double sparsity": weak temporal variation and sparse spatial distribution, thus more representative to reveal long-term latent risks.

[0041] In order to avoid setting too rough amplitude limit manually, the present application constructs a disturbance amplitude threshold A0 with self-adaptability based on statistical method and data distribution characteristics. The construction method is as follows: For all periodic feature points ΔS(t) sequences, the standard deviation σ_RMS of the change rate is calculated as a reference quantity of the overall disturbance level.

[0042] The disturbance threshold A0 is defined as k times σ_RMS, that is, A0= k × σ_RMS. Wherein k is an adjustment factor, the recommended value is between 0.3 and 0.7, which is adjusted according to the application scene, the lower the value of k, the smaller the selected disturbance.

[0043] After calculating the mean and covariance of the one-dimensional feature set composed of all RMS values, the Mahalanobis distance formula is used to judge the abnormal disturbance points. The points with Mahalanobis distance greater than the set threshold M0 (such as 2.5) will be removed to avoid high-amplitude disturbance points from entering the D1 set.

[0044] Verify whether the disturbance distribution of all points in D1 obeys the normal distribution, if not, perform Box-Cox transformation on the RMS value to ensure the uniformity of the disturbance scale within a controllable range. ​

[0045] 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.

[0046] 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.

[0047] 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: 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.

[0048] 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.

[0049] 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.

[0050] 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.

[0051] 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: 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.

[0052] The autocorrelation function method is used to analyze the periodicity of S*(t), and the time lag τ corresponding to the most significant autocorrelation peak is identified 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 temperature and water level.

[0053] By finding the time interval between the first significant change point t0 of the curve and its maximum response point t_max, defined as the delay time Δt = t_max − t0. The delay time can be used to evaluate the response lag of the engineering structure to external changes. If Δt is long, it may represent abnormal stress conduction.

[0054] Statistical analysis of all maxima and minima in the normalized curve records the number of extreme points N_e, the average distance between extreme points Δx_avg, and the average amplitude difference A_diff_avg, reflecting the instability of the structure response.

[0055] Each data point will correspond to a set of five-dimensional feature vectors, including a total of five scalar values in the above four types of features, which are used to construct its hidden danger behavior pattern.

[0056] After summarizing the feature vectors of all data points in D1, the original hidden danger parameter set P_raw is formed. To ensure the consistency of subsequent model input, the set needs to be standardized and dimensionally reduced, with the following specific steps: For each feature dimension in P_raw, the Z-score method is used for standardization, and the standardization formula is: Z = (X− μ) ÷ σ; where X is the original feature value, μ is the mean of this dimension, and σ is the standard deviation. After standardization, all features have a standard distribution with zero mean and unit variance.

[0057] The processed parameter set is named as the hidden danger feature parameter set P, which is used for subsequent coupling mapping with the spatial model. The parameter set is a matrix structure, with the number of rows equal to the number of data points in D1 and the number of columns equal to the number of feature dimensions after standardization.

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

[0059] When collecting monitoring data, each sensor or image acquisition point has a clear spatial coordinate (x, y, z), which is derived from the engineering layout drawing or BIM model. When forming the feature parameter set P, it needs to be one-to-one corresponding to the corresponding coordinates, constructing a feature point set Ps with spatial labels, where each element includes: Point number, spatial position (x, y, z), change slope, fluctuation period, response delay time, local extreme value number, extreme value amplitude difference, and other characteristic parameters.

[0060] Since the digital twin model is usually established on an engineering coordinate system (such as geographic coordinates or structural design coordinates), there may be offsets or precision differences in the collected coordinates, so coordinate unification is required. Preferably, a rigid transformation matrix or ICP (Iterative Closest Point) algorithm is used to align the original coordinate system to the twin model coordinate system.

[0061] After completing the coordinate mapping, each feature point is written into the database as an attribute item of the twin model grid (or structural unit). If the twin model is based on structural grids as the basic unit, the point data can be attributed to the corresponding unit through a spatial indexing algorithm (such as KD tree). Finally, each grid element in the model can contain the following information: Grid ID, number of feature points, mean of main feature value, maximum slope, minimum period, risk weight value, etc.

[0062] According to the main feature parameters (such as change slope or delay time), the point data or grid data is rendered in the three-dimensional model using color gradient, forming a preliminary spatial thermal distribution map. The larger the risk index value, the darker the color, providing a visual basis for subsequent regional aggregation analysis.

[0063] After the feature points are injected into space, further clustering analysis of their distribution and feature consistency is required to find local areas with trend aggregation. This process mainly includes the following steps.

[0064] For each feature point mapped into the model, a fixed radius sphere neighborhood method is used to define its neighborhood radius R (for example, 5 meters). Search for other feature points within its sphere to form a neighborhood set. This operation can be accelerated using KD tree or octree.

[0065] For each feature point and its neighborhood data, the following two indicators are calculated: Neighborhood high-risk point density p: defined as the proportion of points that meet a certain feature (such as change slope greater than threshold k0); Neighborhood feature consistency s: defined as the standard deviation of a certain main feature (such as fluctuation period) of all points in the neighborhood.

[0066] If the density p in a neighborhood is higher than the threshold p0 (such as 0.7) and the standard deviation s is lower than s0 (such as 0.1), then the neighborhood is considered as an abnormal region candidate with "feature trend consistency".

[0067] All points meeting the above conditions are input into the clustering algorithm for spatial aggregation. Preferably, the density-based clustering method DBSCAN is used, with the minimum number of points minPts (e.g., 5) and the maximum radius ε (e.g., 3 meters) set, to extract dense and highly consistent regions in space.

[0068] Each clustering result represents a possible hidden hazard trend aggregation area, is given a unique identifier ID, and extracts its boundary point coordinates, the number of included feature points, and the main feature statistical mean value.

[0069] The final identified aggregation area needs to be clearly marked in the digital twin model and used as input for subsequent simulation calculation and inspection instruction scheduling.

[0070] Each clustering area is given a unique number (e.g., R_001, R_002, etc.) and the following information: region voxel boundary coordinates, number of included feature points, maximum risk indicator value, and average delay time.

[0071] A "potential abnormal area identifier" field is added to the structure unit data structure of the twin model, and the region number is bound to the covered unit, forming a spatial-attribute-risk multidimensional mapping structure.

[0072] The complete attribute information of each abnormal area is saved in the model database or nested structure body.

[0073] The marked area is presented in the twin model as a semi-transparent red area or dynamic heat color block through a three-dimensional visualization tool (such as Unity or Cesium), enabling engineers to visually identify and interact with the risk area.

[0074] In each potential abnormal area identified, there are multiple feature points, each with historical monitoring time series data. To reflect the overall structural response evolution trend of the area, curve modeling of these data is needed.

[0075] First, extract the time series data of all feature points in the potential abnormal area from the database, including at least: Structural strain time series: reflects the deformation of the structure under external load or internal force; Pore water pressure time series: reflects the internal seepage or pressure change state of the dam body; Microseismic response time series: contains vibration amplitude and frequency information, reflecting subtle structural disturbances; Temperature change time series: used to eliminate false responses caused by thermal expansion and contraction; Inspection image analysis indicator sequence: such as the change of crack length and width over time.

[0076] To improve the accuracy of trend identification, the following processing is required for the original time series: Polynomial trend line or moving average method is used to remove long-term linear drift; Time synchronization is performed on data with different sampling frequencies to unify the sampling interval (e.g., once per hour); Missing data is filled using cubic spline interpolation, and outliers are replaced with median values; Physical quantities are normalized to ensure dimensional consistency and facilitate feature comparison.

[0077] The processed data is applied to local weighted regression (LOWESS) for smoothing to obtain continuous response curves. Each physical quantity forms a structural change trend curve, which together forms a trend curve group for the region.

[0078] An association matrix is constructed between the trend curves, such as calculating the Pearson correlation coefficient between the strain curve and the pore pressure curve, and the mutual information index between the microseismic frequency and crack propagation, etc., to identify which physical changes are coupled and driven, and to determine whether the damage is caused by single factor or multiple factors.

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

[0080] 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; seasonal ground temperature change trend. All boundary data are aligned to the trend curve time axis according to a uniform time step.

[0081] The response coupling relationship between the trend curve and the boundary condition is analyzed by the following indicators: Pearson correlation coefficient (r): used to determine whether it is a linear driving relationship; Mutual information value (MI): used to determine the strength of non-linear driving relationship; Time lag analysis value: used to identify lag response, such as strain delay caused by temperature change.

[0082] If the structural response curve and the specific boundary condition show significant correlation (such as r greater than 0.6 or MI value greater than 0.3), it can be inferred that the response of this region is closely related to this boundary factor.

[0083] According to the existing micro-damage database, load typical micro-damage response templates, including but not limited to: Opening crack mode: characterized by rapid increase in strain and decrease in pore pressure; Shear dislocation mode: characterized by microseismic frequency increase and strain change in multiple directions simultaneously; Seepage scouring mode: characterized by frequent pore pressure fluctuations and enhanced microseismic amplitude; Structural relaxation mode: manifested as response hysteresis, low-frequency amplitude rise.

[0084] By the least mean square error or maximum correlation coefficient method, the trend curve group is matched with the typical template, and the micro-damage type that best fits the region is obtained.

[0085] According to 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 structural strain tensor direction, the defect action direction and evolution direction are judged to form a spatial description attribute.

[0086] The estimated micro-damage type and its spatial position, action direction, etc. information are structured and arranged to generate the inversion defect feature set Q.

[0087] Each defect unit contains the following fields: Defect type (opening, shear, scour, etc.); Spatial position (x, y, z); Defect main direction (such as X-axis negative direction, 45° direction of slope, etc.); Damage intensity index (such as maximum strain increment, microseismic energy, etc.); Damage trend description (rise, stable, sharp fluctuation, etc.).

[0088] All defect information is saved in a structured list or JSON object form.

[0089] All defect points form the inversion defect feature set Q, which is an important input for the initial defect boundary conditions in the subsequent local simulation model, and also serves as incremental information for the dynamic evolution and update of the twin model.

[0090] The inversion defect feature set Q is obtained from the previous steps, recording the position coordinates, type, action direction, risk intensity, and trend parameters of each potential micro-damage.

[0091] For each defect record, its spatial position (x, y, z) is taken as a geometric embedding point, combined with its damage type, to construct the corresponding structural disturbance entity. For example: Crack-type defects are modeled as "zero-thickness interfaces" and realized through the definition of contact surface elements; Material softening-type defect regions are disturbed by reducing the elastic modulus or introducing damage constitutive relations; Pore pressure anomaly-type defects are realized by setting local high-pressure sources or permeability sudden change zones in the seepage field.

[0092] A local simulation sub-area is constructed around each defect point, and the spatial range is usually 3 to 5 times the defect influence radius. The mesh division preferably adopts tetrahedral or hexahedral finite elements, and anisotropic mesh encryption is performed according to the main direction of the defect to enhance the simulation accuracy of the stress concentration area.

[0093] The current measured temperature field, water pressure field and structural load are input as initial fields, and the boundary conditions in the future preset period (such as 7 days, 30 days or 90 days) are set, including: external water level change curve; environmental temperature prediction curve; load change simulation (such as upstream water storage, flood discharge, seasonal deformation, etc.).

[0094] The boundary condition data can be provided by historical model prediction, engineering operation plan or weather forecast model.

[0095] On the constructed local model, a complete set of multi-field coupled control equations is established, and time step solving is performed.

[0096] The present application adopts the following three groups of control equations: Heat conduction equation: the control equation of temperature T changing with time t is: specific heat multiplied by density multiplied by the first order derivative of temperature with respect to time, equal to the thermal conductivity coefficient multiplied by the second order derivative of temperature with respect to spatial coordinates.

[0097] Seepage control equation: pore pressure P changes in accordance with Darcy's law and continuity equation, which can be expressed as: pore pressure changes with time equal to the product of permeability coefficient and pressure gradient, and is affected by structural deformation.

[0098] Stress balance equation: the structure stress σ follows the force balance condition, that is: the derivative of internal stress tensor with respect to spatial coordinates is equal to the body force density, combined with the constitutive relation (such as Drucker-Prager model) to reflect the nonlinear response of materials.

[0099] The above three fields interact through coupling terms, for example: temperature affects material expansion, thereby causing stress change; seepage field affects effective stress of structure; stress deformation in turn affects porosity, and then affects seepage path. The coupling strategy can be completed by "step-by-step iterative coupling" or "strong coupling element method". The present application preferably adopts the element strong coupling method to reduce error transmission.

[0100] The finite element method is adopted for spatial discretization, and the backward Euler format is adopted for time dimension. Set the time step Δt (such as 1 hour), the total simulation time T (such as 720 hours), a total of 720 steps are required. After each step solving, the principal value of stress tensor of each structure element in the region is recorded, and time series data is formed.

[0101] After the simulation is completed, the stress change curve with time of the key position in the potential abnormal area, i.e., the stress response evolution curve, is extracted.

[0102] The original output may have discrete oscillation fluctuations, which need to be processed by the following method: The sliding window average method (window width is 5-10 time steps) is used; First-order low-pass filtering is performed to retain the main trend; Finally, the continuous and readable stress evolution curve S(t) is obtained.

[0103] The following indexes are extracted from S(t): Maximum principal stress σ_max; Average growth rate (stress growth per unit time); Stress growth acceleration (i.e., the slope of the slope); Stress trend fitting type (linear, power function, exponential growth, etc.); The trend type can be optimized according to the fitting residual, and the commonly used fitting function is: ; Where a, b, and c are fitting parameters, and b reflects the acceleration degree.

[0104] The curve is embedded in the digital twin system interface in the form of a two-dimensional chart, and the trend acceleration area is marked with a risk color label on the three-dimensional model.

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

[0106] After the local multi-field simulation is completed, the stress principal value time series curve S(t) of the key point in the potential abnormal area is obtained, which is expressed as: , wherein is the maximum principal stress, t is the time, and n is the total number of time steps. To comprehensively evaluate the change trend and stability of the curve, the present application extracts key features from the following three dimensions: The first derivative of the curve S(t), i.e., the stress change rate per unit time, is calculated, and the formula is: R=maximum stress growth rate ÷ average stress change rate; wherein the maximum growth rate is the maximum value of 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 in an amplification state.

[0107] The mean difference between adjacent windows is calculated using the sliding window method (window length W, generally set to 5-10 time steps), and if a difference Δσ is greater than a set threshold value (such as 2 times the reference fluctuation amplitude), it is identified as a "step event".

[0108] After accumulating all the step events, define the step intensity index S as: S = total number of step events × average step amplitude; this index reflects whether there is a sudden change in structural response, which may be caused by boundary disturbance, structural damage or critical instability.

[0109] Perform a fast Fourier transform (FFT) or continuous wavelet transform on S(t) to obtain its frequency component spectrum. After identifying the frequency mutation points (such as main frequency jump, frequency band expansion phenomenon), define the frequency disturbance factor as: F_var = frequency mutation segment energy proportion × main frequency offset amplitude; this factor is used to reflect whether there is an abnormal disturbance in the frequency domain in the structural response, which is usually associated with physical mechanisms such as material softening and boundary degradation.

[0110] After obtaining the three basic indexes R, S, and F_var, the invention uses a weighted linear combination to construct a risk score value V: ; wherein: is a weight coefficient, satisfying ; the weight can be set according to engineering experience or optimized through historical case training, and the example setting is ; V is the final risk score value, which is specific according to the normalization method.

[0111] By analyzing typical engineering accident cases and daily operation data, the corresponding V value distribution is counted, and the scoring space is automatically divided using K-means clustering or piecewise regression method. Finally, the following threshold values are set: T1: the boundary between low risk and medium risk, recommended to be 0.25; T2: the boundary between medium risk and high risk, recommended to be 0.50; T3: the boundary between high risk and extremely high risk, recommended to be 0.75.

[0112] Users can also adjust according to actual project experience to realize the migration adaptability of the model in different projects.

[0113] When the calculated risk score value V satisfies the following conditions, output the corresponding risk level L: If V ≤ T1, it is "low risk"; If T1 < V ≤ T2, it is "medium risk"; If T2 < V ≤ T3, it is "high risk"; if V > T3, it is "extremely high risk".

[0114] Each risk level will be accompanied by trend information as a label, including: rising trend (R continuously increasing); fluctuation trend (S high frequency); frequency disturbance enhancement (F_var sudden increase); Multiple factors superposition (R, S, F_var are high).

[0115] Risk level is mapped to the surface of the digital twin three-dimensional model in a color gradient manner. For example: Green represents low risk; Yellow represents medium risk; Orange represents high risk; Red or purple represents extremely high risk.

[0116] The color mapping follows a linear or exponential scaling relationship, making it easy to visually determine the structure state.

[0117] Trend labels (such as "rapid growth", "frequency disorder", "multiple factor superposition") are displayed on the surface of the corresponding area or in a floating window in the form of arrows, labels or dynamic icons.

[0118] Each risk score value V, level L and trend description is written into the structure unit attribute table of the twin model, which is used for model state version management and evolution path comparison, and can trigger the following operations: Risk alert push; automatically generate inspection or detection suggestions; start emergency response or simulation prediction tasks.

[0119] Example 2, please refer to Figure 2 The water conservancy project supervision analysis system based on digital twin described in the embodiment includes: Data acquisition module: acquire multi-source monitoring data of target water conservancy project in operation period, and synchronize and align the data according to spatial position and time sequence to construct standardized monitoring data set D; Data subset establishment module: based on the historical change law in data set D, identify data points with periodic fluctuation characteristics, and extract a data subset D1 with weak disturbance characteristics combined with 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 the local area with abnormal evolution trend, and mark it as a potential abnormal area in the digital twin model; Inversion defect module: based on the corresponding historical time series data in the potential abnormal area, construct the local structure change trend curve of the potential abnormal area, and inversely deduce the possible types and distribution directions of micro-damage combined with the current environmental boundary conditions to obtain the inversion defect feature set; The coupling calculation module inputs the inversion defect feature set into the digital twin model, performs local thermal-infiltration-force coupling calculation, and obtains a stress response evolution curve of the potential abnormal area in a future preset period; The risk level division module calculates a risk score value V of the potential abnormal area according to a growth rate, a step feature and a frequency mutation in the stress response curve, and outputs a risk level and corresponding evolution trend information of the potential abnormal area if the V exceeds a set threshold T.

[0120] The above merely provides a specific implementation of the present application, but the protection scope of the present application is not limited thereto, and any person skilled in the art can easily think of changes or replacements within the technical range disclosed by the present application, which should be covered within the protection scope of the present 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; 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 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, a local structure change trend curve of the potential abnormal area is constructed, and the possible types and distribution directions of micro-damage are inversely deduced combined with the current environmental boundary conditions to obtain an inversion defect characteristic set; The inversion defect characteristic set is input into the digital twin model, local thermal-infiltration-force 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 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 surface temperature, structural strain, pore water pressure, microseismic amplitude and frequency, and crack features in artificial inspection images.

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

4. The hydraulic engineering supervision analysis method based on digital twinning according to claim 3, characterized in that: The multi-parameter joint fitting of the data subset D1 comprises: For each data point in D1, its complete time series is extracted, and the missing data is completed by 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 regression fitting, the sliding window method is used to further optimize the curve shape, extract the local change trend, and construct the 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, which comprises: The first derivative of the response curve in the whole period is calculated to obtain the change rate function, so as to extract the maximum change slope value and its 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, the autocorrelation function method is used to identify the main period, and the period length is recorded as the fluctuation period characteristic; The time difference between the initial disturbance point and the maximum response point is calculated by the time delay analysis method, and the delay time parameter is extracted; The positions and amplitudes of all local maximum and minimum points in the curve are identified, the number of extreme points, the average spacing, and the amplitude difference are counted, and a local extreme value feature subset is constructed.

6. The hydraulic engineering supervision and analysis method based on digital twinning according to claim 5, characterized in that: The local area with abnormal evolution trend is screened out, including: Bind the spatial position coordinate information (x, y, z) of each hidden danger feature vector corresponding to the data collection stage 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 define 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. If the high-risk point density in a neighborhood exceeds a certain threshold ρ0, and the standard deviation of the eigenvalues in the region is lower than the threshold σ0, it is judged as an abnormal region with consistent trend; A clustering algorithm is used to spatially cluster the point set that meets the conditions, and multiple local risk regions with consistent evolution are extracted; 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 6, characterized in that: Based on the corresponding historical time series data in the potential abnormal region, the local structural change trend curve of the potential abnormal region is constructed, including: Extract the time series data of all feature points in the potential abnormal region within the historical operation period, including strain, pore water pressure, microseismic response amplitude and frequency parameters; After detrending the time series of each type of physical quantity to eliminate long-term steady-state drift, a moving average and local weighted regression method is used to construct a smooth response curve; Through curve fitting, the time variation gradient, variation rate and key turning points of each type of physical quantity are extracted to construct a set of multivariate structural response trend curves; Align the multiple physical field response curves along the time axis and establish a correlation matrix to identify the cooperative relationship between the changes of different physical quantities.

8. The hydraulic engineering supervision analysis method based on digital twinning according to claim 7, characterized in that: The possible types and distribution directions of micro-damage are inferred by combining the current environmental boundary conditions, including: Obtain the current environmental boundary conditions corresponding to the abnormal region, including air temperature, water level, reservoir water pressure change rate, rainfall, and seasonal ground temperature gradient; Perform correlation analysis on the local structural change trend curve and the current boundary conditions, and use the Pearson correlation coefficient and mutual information to quantify the influence intensity; Combine the multiple typical response modes in the preset micro-damage mechanism library to perform pattern matching on the current trend curve; Select the most possible micro-damage type by the least residual matching principle, and record its occurrence direction and spatial offset position.

9. The hydraulic engineering supervision analysis method based on digital twinning according to claim 8, characterized in that: Perform local thermal-hydraulic-mechanical coupling calculation, including: Set the spatial boundary range and time scale of the coupling solution region, including the defect point and its surrounding structure units with an influence radius of not less than 3 times; Import the temperature, water pressure, and structure load at the current time as the initial boundary conditions, and construct the time step driving boundary by combining the future preset period of meteorological and operation prediction data; 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; The principal value change of the stress tensor of the local defect area is extracted at each time step to form a stress response evolution curve, which is used to represent the structure evolution trend.

10. A water conservancy project supervision and analysis system based on digital twinning, used to implement the water conservancy project supervision and analysis method based on digital twinning in any one of claims 1-9, characterized in that: It comprises: A data acquisition module: collects multi-source monitoring data of the target water conservancy project during the operation period, synchronously aligns them according to spatial position and time sequence, and constructs a standardized monitoring data set D; A data subset establishment module: based on the historical change law in the data set D, the data points with periodic fluctuation characteristics are identified, and combined with their spatial distribution density, a data subset D1 with weak disturbance characteristics is extracted; A 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; A potential abnormal area marking module: the hidden danger feature parameter set is mapped into the spatial coordinate system of the digital twin model, the local area with abnormal evolution trend is screened out, and is marked as a potential abnormal area in the digital twin model; An inversion defect module: based on the corresponding historical time series data in the potential abnormal area, a local structure change trend curve of the potential abnormal area is constructed, and the possible types and distribution directions of micro-damage are back calculated combined with the current environmental boundary conditions to obtain an inversion defect feature set; A coupling calculation module: the inversion defect feature set is input into the digital twin model to perform local heat-seepage-force coupling calculation, and the stress response evolution curve of the potential abnormal area in the future preset period is obtained; A risk level division module: 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 corresponding evolution trend information of the potential abnormal area are output.

Images