Reservoir operation state prediction analysis method based on environmental data

Abstract

The present application relates to reservoir safety monitoring and early warning technical field, especially to a kind of reservoir operating state prediction analysis method based on environmental data.The method comprises the following steps: through the vertical resistivity probe array of layout, using low-frequency alternating current signal acquisition reservoir bottom mud different depth original resistance response sequence;Original resistance response sequence is subjected to polarization effect compensation processing, and the bottom mud resistivity layer spectrum diagram is obtained;Reservoir environmental data are collected;Reservoir environmental data and bottom mud resistivity layer spectrum diagram are subjected to time series coupling analysis, and the bottom mud stability influence factor spectrum is obtained;According to bottom mud resistivity layer spectrum diagram and bottom mud stability influence factor spectrum, bottom mud layered structure model is constructed;The initial stress distribution map is obtained by calculating the basic stress field of bottom mud layered structure model.The present application is monitored from bottom mud inside to the whole chain technical route of multi-factor comprehensive analysis, and the early warning timeliness, accuracy and reliability of reservoir bottom mud disaster are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of reservoir safety monitoring and early warning technology, and in particular to a method for predicting and analyzing the operational status of reservoirs based on environmental data. Background Technology

[0002] Traditional monitoring systems generally suffer from early warning lag. Most monitoring methods can only capture surface displacement or deformation of sediment, resulting in a time lag between changes in the internal stress state of the sediment and the appearance of obvious surface displacement. This leads to an excessively short early warning window, making it difficult to provide sufficient preparation time for emergency response. Traditional technologies have limited monitoring capabilities for the internal structure of sediment, relying mainly on indirect methods such as surface displacement monitoring and tilt monitoring. This results in blind spots in monitoring the physical state and structural changes of deep sediment, making it impossible to grasp early signs of sediment instability in real time. Existing methods typically analyze environmental factors in isolation, ignoring the interactions between multiple environmental factors and their comprehensive impact on sediment stability. This leads to insufficient accuracy in early warning judgments under complex environmental conditions, making it difficult to establish a precise correlation model between environmental changes and sediment stability.

[0003] In summary, existing technologies suffer from problems such as insufficient early warning timeliness, limited internal monitoring capabilities, and inadequate comprehensive analysis of environmental factors, which urgently need to be addressed. Summary of the Invention

[0004] Therefore, it is necessary to provide a method for predicting and analyzing the operational status of reservoirs based on environmental data in order to solve at least one of the aforementioned technical problems.

[0005] To achieve the above objectives, a method for predicting and analyzing the operational status of reservoirs based on environmental data includes the following steps:

[0006] Step S1: Using a vertical resistivity probe array, low-frequency AC signals are used to collect the original resistivity response sequences of the reservoir bottom sediment at different depths; polarization effect compensation processing is performed on the original resistivity response sequences to obtain the bottom sediment resistivity spectral map.

[0007] Step S2: Collect reservoir environmental data; perform time-series coupling analysis on the reservoir environmental data and sediment resistivity spectral diagram to obtain the sediment stability influencing factor spectrum;

[0008] Step S3: Construct a sediment layered structure model based on the sediment resistivity stratification diagram and sediment stability influencing factor spectrum; perform basic stress field calculation on the sediment layered structure model to obtain an initial stress distribution map; perform pore water pressure calculation based on the initial stress distribution map to obtain an effective stress distribution map; calculate the sediment strength ratio distribution map based on the effective stress distribution map; perform stress gradient anomaly calculation based on the effective stress distribution map to obtain a stress gradient anomaly region map; perform potential sliding surface stress evolution analysis on the stress gradient anomaly region map and the strength ratio distribution map to obtain the sediment interlayer stress field spectrum.

[0009] Step S4: Based on the stress field spectrum of the bottom sediment layers, assess the probability of reservoir instability to obtain a regional instability risk map; based on the regional instability risk map, conduct a catastrophic criticality early warning to obtain a reservoir bottom sediment disaster early warning report.

[0010] This invention improves the quality of raw data through high-precision deployment and low-frequency measurement. Precise polarization and temperature compensation eliminate interference factors. High-resolution depth and spatial dynamic interpolation construct a resistivity tomographic map reflecting the internal structure of the sediment, enabling precise, dynamic, and non-destructive monitoring of the sediment's physical state. Rigorous preprocessing of multi-source environmental data ensures data quality. Lag correlation analysis reveals the transmission lag of environmental impacts. Partial correlation and interaction pattern analysis delve into the complex interactions between environmental factors. Critical conditions are extracted by combining historical events, establishing a more accurate quantitative correlation model between environmental changes and sediment stability, enhancing the scientific rigor of predictions. Based on resistivity inversion of physical parameters and identification of weak layers, a three-dimensional model accurately reflecting the heterogeneous structure of the sediment is constructed. A fine finite element mesh and transient seepage model accurately calculate the internal stress field and pore water pressure, especially excess pore pressure. Multi-dimensional analysis and time evolution tracking of stress gradients enable the identification of micro-stress anomalies before macroscopic displacement. Combined with dynamic assessment of the potential sliding surface safety factor, a comprehensive and dynamic monitoring map of the sediment's internal mechanical state is provided, significantly advancing the detection time of instability signs. Establishing a historical database of precursors to disasters provides a benchmark for current anomaly detection. Assessments based on multiple models and historical similarity improve the accuracy of instability probability judgments. Considering spatial correlation and cascading effects, the risk of disaster spread is assessed, and the time window for risk outbreaks is predicted. Finally, by calculating graded critical indices and generating structured early warning reports, complex analysis results are transformed into intuitive and timely early warning information, providing early, comprehensive, and scientific decision-making basis for reservoir safety management and effectively improving reservoir operational safety. Therefore, this invention provides a reservoir operational status prediction and analysis method based on environmental data. By constructing a sediment internal state monitoring architecture based on resistivity tomography, real-time monitoring of sediment microstructure changes is achieved. The finite element method is used to reconstruct the internal stress field distribution of the sediment, capturing early stress anomaly signals. A time-series coupling model of multi-source environmental factors and sediment stability is established to comprehensively analyze the impact of environmental changes on sediment stability. This full-chain technical route, from sediment internal monitoring to multi-factor comprehensive analysis, significantly improves the timeliness, accuracy, and reliability of early warnings for reservoir sediment disasters. Attached Figure Description

[0011] Figure 1 This is a flowchart illustrating the steps of a method for predicting and analyzing the operational status of a reservoir based on environmental data.

[0012] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0013] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0014] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0015] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0016] In this embodiment of the invention, reference Figure 1 The diagram shown is a flowchart illustrating the steps of the reservoir operation status prediction and analysis method based on environmental data according to the present invention. In this example, the reservoir operation status prediction and analysis method based on environmental data includes the following steps:

[0017] Step S1: Using a vertical resistivity probe array, low-frequency AC signals are used to collect the original resistivity response sequences of the reservoir bottom sediment at different depths; polarization effect compensation processing is performed on the original resistivity response sequences to obtain the bottom sediment resistivity spectral map.

[0018] In this embodiment of the invention, key monitoring areas of the reservoir are determined using high-precision positioning equipment and historical sedimentary data. An underwater robot is used to deploy a vertical resistivity probe array containing 10 to 15 measuring electrodes (depth intervals of 10 to 20 cm). The probes are calibrated using standard resistance blocks to obtain a probe position depth table containing spatial coordinates and electrode depths. Based on the probe position depth table, a multi-channel measuring instrument sends low-frequency AC signals of 0.1 to 10 Hz and 5 to 20 mA to the probe array. A four-electrode method is used to repeatedly measure each depth point (5 times, 100 Hz sampling, lasting 30 seconds) to obtain the raw resistance response sequence containing voltage / current time series. The raw resistance response sequence is then subjected to time-domain differential or complex resistivity processing to eliminate electrode polarization effects, and temperature compensation based on empirical formulas is performed using temperature sensor data. A set of corrected resistivity values ​​was obtained; cubic spline interpolation was used to generate a vertical depth distribution curve with a depth resolution of 1 cm from the set of corrected resistivity values; finally, three-dimensional kriging space interpolation (horizontal resolution 5 m × 5 m) was performed on the vertical curves at different probe positions, and time series data (1 hour resolution) were integrated to construct a sediment resistivity stratum map containing three-dimensional information of time, spatial location and depth.

[0019] Step S2: Collect reservoir environmental data; perform time-series coupling analysis on the reservoir environmental data and sediment resistivity spectral diagram to obtain the sediment stability influencing factor spectrum;

[0020] In this embodiment of the invention, multi-source environmental data, including reservoir water quality (pH, dissolved oxygen, turbidity, water temperature), meteorological data (rainfall, air temperature, air pressure), and hydrological data (water level, inflow, outflow), are collected. Time alignment (hourly intervals), missing values ​​are filled with nearest neighbor means, outliers are removed using median filtering, and Z-score standardization is performed to obtain a time series table of multi-source environmental parameters. The lag correlation coefficients (0 to 72-hour lags) between the time series of each environmental parameter and the resistivity time series of each depth layer (1 cm resolution) of the sediment are calculated. The correlation stability is assessed using a 7-day sliding window. STL decomposition of seasonal parameters is performed, and residual correlation is analyzed to obtain a factor-depth response matrix containing the degree of influence, lag time, and stability score. Based on... Factor-depth response matrices are used to construct factor correlation networks. Partial correlation analysis is employed to identify synergistic or antagonistic direct interaction patterns among environmental factors. The time lag of the impact of these interaction patterns on resistivity is analyzed. Critical conditions or thresholds for key interaction patterns are determined by combining historical instability events. The performance differences of these interaction patterns at different depth layers are analyzed. Finally, a factor interaction effect map is drawn to intuitively show the factor interactions, critical conditions, and depth distribution. Finally, by combining the factor interaction effect map, factor-depth response matrix, and historical instability event data, decision trees or logistic regression models are used to extract the influence weights of each environmental factor (including single factors and combined factors) on sediment stability and the critical threshold ranges that trigger instability, forming a sediment stability influencing factor spectrum.

[0021] Step S3: Construct a sediment layered structure model based on the sediment resistivity stratification diagram and sediment stability influencing factor spectrum; perform basic stress field calculation on the sediment layered structure model to obtain an initial stress distribution map; perform pore water pressure calculation based on the initial stress distribution map to obtain an effective stress distribution map; calculate the sediment strength ratio distribution map based on the effective stress distribution map; perform stress gradient anomaly calculation based on the effective stress distribution map to obtain a stress gradient anomaly region map; perform potential sliding surface stress evolution analysis on the stress gradient anomaly region map and the strength ratio distribution map to obtain the sediment interlayer stress field spectrum.

[0022] In this embodiment of the invention, based on the sediment resistivity spectral diagram and the sediment stability influencing factor spectrum, the improved Archie formula (ρ=a×ρ_w×φ) is used. -m Invert the sediment porosity, water content, and density, and estimate the cohesion c and internal friction angle φ by combining the sediment type and water content. iBy referencing the influence factor spectrum and correcting the intensity parameters, a sediment physical property profile containing the vertical distribution of various physical parameters is obtained. The resistivity profile is differentiated to detect layer peaks and identify candidate points at layer interfaces. Physical property feature vectors are extracted from the candidate points, and K-means clustering (K=3) is performed to stratify the layers, identifying interlayer transition zones and key weak layers with low intensity and high sensitivity. The underwater topographic DEM and physical property stratification scheme are integrated, and a three-dimensional geometric model of each layer interface is constructed through spatial interpolation. Three-dimensional kriging space interpolation is performed on the discrete physical parameters to obtain the stratified physical parameter field. Geometric, parameter, and weak layer information are integrated to generate a sediment stratification structure model. Based on the stratified physical parameters, the parameters (E, ν, c, φ) of the elastoplastic constitutive model are determined. i The layered model is discretized into finite element meshes with horizontal dimensions of 1-5 meters and vertical dimensions of 5-20 centimeters (with weak layers refined to 2 centimeters). Material properties are assigned, and the initial total stress distribution under self-weight and water pressure is calculated. The pore water pressure field is calculated using a dynamic permeability coefficient field (corrected based on resistivity and environmental factors) and a transient seepage model. The effective stress distribution is obtained by subtracting the pore water pressure from the total stress. The effective stress field is analyzed by calculating multi-scale spatial gradient, shear stress curl, and the ratio of stress gradient to local strength, and its temporal evolution trend is analyzed to identify stress gradient anomaly regions. The location and shape of potential sliding surfaces are identified by combining the stress gradient anomaly regions and the shear stress to shear strength ratio map. The safety factor of each potential sliding surface is calculated using the limit equilibrium method (such as the simplified Bishop method). Finally, the real-time state, rate of change, and correlation with environmental factors of effective stress, pore water pressure, potential sliding surfaces, and safety factors are comprehensively analyzed, and future development trends are predicted to generate a sediment interlayer stress field map containing all dynamic mechanical state information.

[0023] Step S4: Based on the stress field spectrum of the bottom sediment layers, assess the probability of reservoir instability to obtain a regional instability risk map; based on the regional instability risk map, conduct a catastrophic criticality early warning to obtain a reservoir bottom sediment disaster early warning report.

[0024] In this embodiment of the invention, historical sediment disaster data (landslides, liquefaction, etc.) of reservoirs are collected, and environmental and monitoring data before the disasters are retrospectively analyzed. The stress field characteristics before the disasters are simulated or inverted, and typical stress field precursor patterns (such as high intensity ratio, high gradient, high pore pressure ratio, low Fs and decreasing) are extracted to construct a disaster precursor feature library. Based on the current interlayer stress field map of the sediment, and in comparison with the thresholds and patterns in the disaster precursor feature library, abnormal stress field regions are detected, and the anomaly type, degree, and rate of change are recorded to obtain a table of abnormal stress field regions. Potential instability modes (shear sliding, liquefaction, etc.) are identified based on combinations of abnormal features. The safety factor and rate of change of the potential sliding surface are evaluated for the shear sliding mode, and the pore water pressure ratio and material sensitivity are evaluated for the liquefaction mode. Simultaneously, the similarity between the current anomaly and historical precursor cases is calculated. Combining the pattern evaluation results and historical similarity, a probability model is used to calculate the probability of local instability in each abnormal region. The system employs spatial interpolation to obtain the regional instability probability distribution. Considering the interaction and potential cascading effects of stress / pore pressure between regions, it simulates the possible propagation and impact of disasters, assessing the risk of disaster spread. Integrating the spatial distribution of risk and the rate of change of key parameters (such as the rate of change of the safety factor), it predicts the time window for reaching the critical state, forming a regional instability risk map that includes spatial risk distribution and time windows. Based on the highest probability and high-risk area range of the regional instability risk map, as well as the highest anomaly degree and predicted time window of the stress field anomaly area table, it calculates a comprehensive sediment disaster criticality index and delineates warning levels (normal, caution, warning, danger, emergency). Finally, based on the criticality index and regional instability risk map information, it generates a reservoir sediment disaster early warning report containing the warning level, anomaly location, potential disaster type, impact range, predicted occurrence time window, and recommended countermeasures, and disseminates it through multiple channels.

[0025] Most importantly, the polarization effect compensation treatment is specifically as follows:

[0026] Extract the phase difference feature matrix of the original resistance response sequence;

[0027] Frequency response curves are constructed from the original resistance response sequence to obtain a set of frequency response parameters;

[0028] Collect temperature data around the measurement point; generate a temperature distribution profile based on the temperature data;

[0029] Temperature-corrected resistivity set is obtained by performing temperature-effect correction on the frequency response parameter set based on the temperature distribution profile.

[0030] Based on the phase difference characteristic matrix, the electrode contact impedance of the temperature-corrected resistivity set is eliminated to obtain the corrected resistivity value set.

[0031] In this embodiment of the invention, a raw resistance response sequence containing voltage and current time-series data is received. For each set of repeated measurements at each measurement point (specific probe array and depth) in the sequence, a Fast Fourier Transform (FFT) is applied to extract the complex spectral values ​​of voltage and current at the excitation frequency (e.g., set to 5 Hz). The phase difference φ = phase(V(f0)) - phase(I(f0)) between the voltage and current spectra is calculated, and the average phase difference value of multiple repeated measurements is stored as a phase difference feature matrix, which reflects the imaginary part of the impedance characteristics caused by electrode and sediment polarization. Simultaneously, based on the voltage and current spectrum amplitudes extracted by FFT, the apparent impedance amplitude |Z| = |V(f0)| / |I(f0)| corresponding to the excitation frequency is calculated, and the average amplitude of multiple repeated measurements is stored as a frequency response parameter set. Next, real-time temperature data at different depths around the measurement point is acquired using a temperature sensor integrated on the probe array, and a sediment temperature distribution profile covering the entire measurement depth is generated using an interpolation method (such as linear or cubic spline interpolation). Based on the temperature distribution profile, and according to the temperature T at the time of measurement, an empirical temperature correction formula (e.g., ρ) for saturated sediment is used. 25 =ρ_a / [1+0.02×(T-25)], where ρ_a is the apparent resistivity calculated using the apparent impedance amplitude and geometric factor. Temperature correction is applied to the apparent resistivity in the frequency response parameter set to obtain a temperature-corrected resistivity set, which represents the apparent resistivity at standard temperature. Finally, to eliminate the influence of electrode polarization and sediment induction polarization on the resistivity amplitude, complex resistivity theory is used, based on the phase difference φ in the phase difference characteristic matrix and the temperature-corrected apparent resistivity ρ in the temperature-corrected resistivity set. 25 Calculate the true resistivity value ρ_corrected=ρ 25 ×cos(φ). This operation extracts the real part of the complex impedance, resulting in a set of corrected resistivity values ​​that more closely approximate the bulk resistivity of the sediment.

[0032] Preferably, step S1 includes the following steps:

[0033] Step S11: Obtain and deploy vertical resistivity probe arrays in key monitoring areas of the reservoir based on the coordinates of the reservoir bottom sediment area and historical sedimentation data. Each probe array contains 10-15 measuring electrodes at different depths with a depth interval of 10-20 cm, and obtain a probe position depth table.

[0034] Step S12: Send a low-frequency AC signal of 0.1-10Hz to the probe array according to the probe position depth table. The signal strength is 5-20mA. Then perform 5 repeated measurements at each depth point with a sampling rate of 100Hz and each measurement lasts for 30 seconds. Finally, obtain the original resistance response sequence.

[0035] Step S13: Perform polarization effect compensation on the original resistance response sequence to obtain a set of corrected resistivity values;

[0036] Step S14: Perform depth interpolation on the set of corrected resistivity values ​​to obtain the depth distribution curve;

[0037] Step S15: Construct a resistivity stratum map of the bottom sediment based on the depth distribution curve.

[0038] In this embodiment of the invention, a high-precision differential GPS receiver is first used to accurately locate the reservoir and obtain WGS84 coordinate data of the reservoir's bottom sediment area. Simultaneously, geological exploration reports from the reservoir's construction period, historical sediment deposition measurement reports, and historical landslide accident investigation reports are reviewed to extract the average thickness, layering characteristics, distribution of major soil types, and locations of historically unstable areas of the bottom sediment. By comprehensively analyzing the aforementioned coordinate data and historical sedimentary data, key monitoring areas are identified in the reservoir's bottom sediment area, specifically those near important structures (such as dams and spillways), known historically unstable areas, or areas with high sediment deposition rates. In each key monitoring area, an underwater robot carrying a vertical probe is used to precisely insert a vertical resistivity probe array into the bottom sediment at a specified depth. Each probe array is made of a rigid insulating material (such as a high-strength PVC pipe), on which 12 annular stainless steel measuring electrodes are fixed. The electrodes are evenly spaced along the vertical direction, with the electrode spacing strictly set at 15 cm, thereby achieving layered monitoring of the bottom sediment up to a depth of 1.8 meters below the surface. After the probe array is deployed, the electrode sequence of each probe array is calibrated using a standard resistance block (e.g., a standard tank filled with a KCl solution with a known conductivity of 1413 μS / cm). The probe arrays are sequentially immersed in the standard tank, and the resistance of each electrode pair (e.g., using a four-electrode method) is measured under standard conditions using a precision resistivity meter. The deviation between the measured values ​​and the theoretical resistance values ​​of the standard resistance block is recorded, and a calibration coefficient table is generated for subsequent data correction. Finally, a probe position depth table is formed by recording the unique identifier of each deployed probe array, its precise latitude and longitude coordinates, and the accurate depth of each measuring electrode relative to the sediment surface. A dedicated multi-channel resistivity meter is connected to the deployed vertical resistivity probe array via an underwater cable. The meter automatically controls electrode switching based on the electrode depth information recorded in the probe position depth table, using a four-electrode method for measurement. For example, using a Wenner arrangement, current is transmitted through the outer pair of electrodes (A and B), and voltage is measured through the inner pair of electrodes (M and N). The measuring instrument outputs a sinusoidal AC constant current signal with a frequency of 5Hz and an intensity of 10mA. For different electrode combinations on the probe array (e.g., selecting electrodes 1 and 4 as current transmitting electrodes A and B, and electrodes 2 and 3 as voltage measuring electrodes M and N to measure the resistivity of the corresponding depth layer; then switching electrode combinations to measure other depth layers), the measuring instrument performs 5 independent measurement cycles for each specific electrode combination. In each measurement cycle, the signal transmission lasts for 30 seconds, and the measuring instrument synchronously acquires the time series data of the transmitted current and measured voltage at a sampling rate of 100Hz. The voltage time series acquired in each measurement cycle is divided by the current time series to obtain the time series of the original resistance values. The time series data of the original resistance values ​​of all measurement points and measurement cycles are collected, and the corresponding timestamps, probe array IDs, and electrode combination IDs are recorded to form the original resistance response sequence.The system receives raw resistance response sequence data. Voltage and current time-series data for each measurement point and each measurement cycle are processed to analyze the phase difference between the voltage and current signals. Due to electrode polarization, frequency-dependent impedance components are generated in actual measurements. Using the time-domain difference method, the instantaneous rate of change of the voltage time series relative to the current time series within 30 seconds of each measurement cycle is calculated, and its decay characteristics are analyzed. Alternatively, using the principle of complex resistivity, the virtual resistance component caused by electrode polarization is identified and separated, and the true real resistance component, i.e., the resistance value after deducting the polarization effect, is extracted. Simultaneously, real-time temperature data of the sediment measurement depth is collected using a temperature sensor integrated on the probe array. Based on the empirical formula for temperature-resistivity of water-saturated sediment, temperature compensation is applied to the resistance value after deducting the polarization effect, correcting it to the resistivity value at a standard temperature of 25℃. The temperature compensation formula is: [Formula omitted for brevity]. Where ρ 25 It is the resistivity at 25℃ (unit: ohm-meter). is the resistivity at temperature T (unit: °C) (unit: ohm-meter), and 0.02 is the empirical temperature coefficient (unit: °C⁻¹). The actual resistivity values ​​at each measurement point, after polarization and temperature compensation, are correlated with the corresponding timestamp, probe array ID, and measurement depth to form a set of calibrated resistivity values. The calibrated resistivity value set is received. For each probe array's resistivity data set at a specific time point (this set contains the calibrated resistivity values ​​of 12 electrodes at different depths on that probe array), cubic spline interpolation is used to interpolate the resistivity values ​​at these discrete depth points. Using the sediment surface as the zero depth point, interpolation is performed vertically downwards to generate a continuous resistivity-depth curve from the sediment surface to the deepest electrode depth (e.g., 1.8 meters). The interpolation process ensures that the curve passes through the actual values ​​of all measurement points and transitions smoothly between points. The interpolation results are output as resistivity values ​​at fixed depth intervals of 1 cm. The interpolation results of each probe array at various time points are saved to form a depth distribution curve dataset containing timestamps, probe array IDs, and resistivity-depth sequences at 1-cm intervals at that probe location. Depth distribution curve data of all probe arrays at different time points are received. This data is organized according to the spatial location (latitude and longitude coordinates) of the probe arrays and time. To construct the three-dimensional resistivity distribution of the sediment region, Kriging interpolation is used to spatially interpolate the discrete probe array data. The entire sediment monitoring area is divided into a grid with a horizontal resolution of 5m × 5m. For each grid cell, the resistivity prediction value of each depth layer (at 1-cm intervals) within that grid cell is calculated using the depth distribution curve data of nearby probe arrays combined with the Kriging interpolation algorithm. The three-dimensional resistivity distribution field data calculated at different time points (e.g., hourly) are stored as a time series. Finally, a resistivity stratum map of the sediment was formed. This map records the resistivity distribution of each depth layer (1 cm resolution) of the sediment at different spatial locations (5 m × 5 m grid resolution) and its variation characteristics over time (1 hour resolution), forming a resistivity distribution dataset containing information in three dimensions: time, spatial location, and depth.

[0039] Preferably, step S2 includes:

[0040] Step S21: Perform environmental data preprocessing on the reservoir environmental data to obtain a time series table of multi-source environmental parameters;

[0041] Step S22: Perform hysteresis correlation analysis on the sediment resistivity spectral diagram and the time series table of multi-source environmental parameters to obtain the factor-depth response matrix;

[0042] Step S23: Perform environmental factor interaction analysis on the factor-depth response matrix to obtain the factor interaction effect diagram;

[0043] Step S24: Extract the critical conditions for stability from the factor interaction effect diagram and the factor-depth response matrix to obtain the spectrum of factors affecting sediment stability.

[0044] In this embodiment of the invention, water quality parameters (e.g., pH value, dissolved oxygen concentration, turbidity, and water temperature are collected using multi-parameter water quality sensors), meteorological data (e.g., rainfall, air temperature, and air pressure are collected using automatic weather stations), and hydrological data (e.g., water level is collected using level gauges, and inflow and outflow are collected using flow meters) are collected using existing automated monitoring stations in the reservoir. These sensors and devices record raw data at a preset frequency (e.g., every 15 minutes or every 30 minutes). After receiving this raw data, time alignment is first performed. If the raw data collection frequency is higher than 1 hour, the data is aggregated into one value per hour using a mean aggregation method; if the raw data collection frequency is lower than 1 hour, the data is interpolated into one value per hour using a linear interpolation method, ensuring that the time series of all environmental parameters are strictly aligned with the time points of the sediment resistivity spectral map, i.e., using hours as the smallest time unit. Next, missing data is filled. For missing data at fewer than 3 consecutive time points, the average value of the values ​​at the nearest adjacent time points is used for filling; for missing data at 3 or more consecutive time points, the data is marked as invalid and not included in subsequent analysis. Outlier removal was then performed. Median filtering with a window length of 5 time points was applied to the time series data of each environmental parameter. Values ​​deviating from the window median by more than 3 times the standard deviation were marked as outliers and removed. Finally, the processed environmental parameter data were standardized using Z-score standardization, converting the time series of each parameter into a distribution with a mean of 0 and a standard deviation of 1. The standardization formula is: Z = (X - μ) / σ, where X is the original value, μ is the mean of the parameter, and σ is the standard deviation. The standardized time series data of each environmental factor were integrated into a table, recording all environmental parameter values ​​for each hour, forming a multi-source environmental parameter time series table. Sediment resistivity tomography (containing resistivity data from different spatial locations, depths, and times) and multi-source environmental parameter time series tables (containing standardized values ​​of different environmental parameters at different times) were obtained. For the sediment resistivity tomography, the resistivity time series of each probe array location (or after spatial averaging) at each depth layer (e.g., every 1 cm depth) was first extracted. Then, the lag correlation coefficients between the time series of each environmental parameter (e.g., rainfall) in the multi-source environmental parameter time series table and the time series of sediment resistivity at each specific depth layer are calculated. The lag time window is set to 0 to 72 hours, i.e., the Pearson correlation coefficients between the environmental parameter time series and the resistivity time series at lags of 0 hours, 1 hour, ..., 72 hours are calculated. The formula for calculating the correlation coefficient is: Where r is the correlation coefficient, xi These are the values ​​of environmental parameters at time point i. It is the mean of the time series of environmental parameters, y i+ τ is the value of resistivity at time point i+τ (lag time). This represents the mean of the resistivity time series, with ∑ indicating summation. For environmental parameters with significant seasonal variations (such as water temperature and water level) or resistivity changes, the STL decomposition method is used to decompose the time series into seasonal, trend, and residual terms. Lagged correlation analysis is then performed on the residual terms to eliminate seasonal interference. A sliding window method (with a window length of 7 days) is used to calculate the stability of the lag correlation coefficient across different time periods; that is, the standard deviation of the correlation coefficient within each window is calculated, with a smaller standard deviation indicating a more stable correlation. Factor-depth pairs with an absolute correlation coefficient greater than 0.5 and a sliding window standard deviation less than 0.1 throughout the monitoring period are selected as combinations with significant and stable correlations. Finally, the maximum correlation coefficient between each environmental factor and each depth layer, the corresponding lag time, and the correlation stability score (e.g., the reciprocal of the standard deviation) are recorded to form a factor-depth response matrix. The rows of this matrix represent environmental factors, and the columns represent sediment depth layers. Based on the factor-depth response matrix, an environmental factor correlation network diagram is constructed. The nodes in the network diagram represent environmental factors, and the lines connecting the nodes represent the correlation strength between factors (which can be the average correlation between factors in the factor-depth response matrix or the correlation at a specific depth). Based on this network diagram, partial correlation analysis is used to calculate the partial correlation coefficient between any two environmental factors, while controlling for the influence of all other environmental factors in the network. Partial correlation analysis can reveal the strength of direct correlations between factors and exclude indirect influences. The calculation results form a factor direct influence matrix. Analyzing the factor direct influence matrix identifies factor combinations with significant positive partial correlations (synergistic enhancement) or negative partial correlations (antagonistic weakening), and classifies these combinations and their interaction types to obtain a factor interaction pattern table. For example, high rainfall and high inflow have a synergistic effect. Based on the factor interaction pattern table and the results of the lag correlation analysis in step S22 (especially the lag time information), the time-lag characteristics of the combined effect of different factor combinations on sediment resistivity are analyzed to obtain a time-lag interaction matrix. This matrix records the key factor combinations and their optimal time lags affecting sediment resistivity. By combining historical monitoring data and known sediment instability events (e.g., from historical disaster records), the numerical characteristics of key factor combinations identified in the factor interaction pattern table before these events were analyzed. Through statistical analysis or machine learning methods (e.g., support vector machine classifiers), the critical threshold ranges or proportional relationships of these factor combinations inducing significant changes in sediment resistivity or instability were determined, resulting in a conditional threshold table. For example, when the ratio of pore water pressure to effective normal stress exceeds a certain threshold, and the rate of water level change exceeds a certain threshold, the risk of instability increases significantly. Using the conditional threshold table and the factor-depth response matrix, the differences in the performance of these interaction patterns and conditional thresholds at different depths were analyzed, identifying the characteristics of specific depth layers being more sensitive to certain factor interaction combinations, resulting in a depth interaction distribution map.Finally, the depth interaction distribution map, factor interaction pattern table, and condition threshold table were integrated to create a factor interaction effect map. This map visually illustrates how different environmental factors interact, their key interaction patterns, critical conditions, and the distribution characteristics of these interactions in the vertical direction of the sediment, thus revealing the complex mechanism by which environmental changes affect sediment stability. The factor interaction effect map (describing the interaction patterns, critical conditions, and depth distribution between factors) and the factor-depth response matrix (describing the influence and lag of individual factors at different depths) were used as inputs. Combined with the time points of sediment instability events marked in historical monitoring data, the characteristics of changes in environmental factors and sediment resistivity tomography in the period preceding these events were retrospectively analyzed. A decision tree learning method or logistic regression model was applied, using standardized environmental factor values ​​(including single factors and combined factors constructed based on the factor interaction pattern table) as input features and the presence of instability signs in the sediment (e.g., rapid decrease in resistivity, decreased safety factor shown in stress field calculations, etc., defined based on historical experience or preset thresholds) as output labels to train a discriminative model. Discriminant rules are extracted from a trained decision tree model. These rules describe the judgment logic for sediment stability when environmental factors reach specific numerical combinations or exceed specific thresholds. Simultaneously, the importance weights of each environmental factor (including single and combined factors) in judging sediment stability are extracted from the model; the weight values ​​reflect the degree of influence of the factor on sediment stability. Based on the discriminant rules and historical data, the critical threshold ranges of each environmental factor are refined. For key interaction combinations in the factor interaction effect diagram, their joint critical conditions (e.g., specific combinations of water level drop rate and pore water pressure change rate) are extracted. Combined with the factor-depth response matrix, the intensity and manifestation of these critical conditions at different depth layers are determined. Finally, a sediment stability influencing factor spectrum is formed. This spectrum records in a structured form the influence weight of each key environmental factor (or factor combination) on sediment stability, its critical threshold range for triggering unstable states (including numerical range or rate of change), and the characteristics or sensitivity differences of the factor (or combination) at different depth layers in the sediment (e.g., shallow, middle, and deep layers). This provides a quantitative influence mechanism model for subsequent sediment internal stress field reconstruction and catastrophic criticality prediction.

[0045] Preferably, the environmental factor interaction analysis in step S2 includes:

[0046] Construct a factor correlation network graph based on the factor-depth response matrix;

[0047] Partial correlation analysis was performed on the factor correlation network diagram to obtain the direct influence matrix of the factors;

[0048] The interaction pattern classification and identification of the direct influence matrix of factors is performed to obtain the factor interaction pattern table;

[0049] Based on the factor interaction pattern table, the time-delay interaction effect is analyzed to obtain the time-delay interaction matrix;

[0050] Determine the conditional threshold table based on the time-delay interaction matrix;

[0051] A deep difference analysis of the factor-depth response matrix was performed based on the conditional threshold table to obtain a deep interaction distribution map;

[0052] A factor interaction effect diagram is drawn based on the deep interaction distribution diagram, the factor interaction pattern table, and the conditional threshold table.

[0053] In this embodiment of the invention, a multi-source environmental parameter time series table is obtained. This table contains standardized values ​​of environmental factors such as water quality parameters (pH, dissolved oxygen, turbidity, water temperature), meteorological data (rainfall, air temperature, air pressure), and hydrological data (water level, inflow, outflow) over a continuous time series. The Pearson correlation coefficient between any two different environmental factor time series in the multi-source environmental parameter time series table is calculated. The correlation coefficient calculation covers the entire monitoring period. An undirected weighted graph is constructed, where nodes represent various environmental factors, and the lines connecting nodes represent the correlation between two factors. The weight of the lines is set to the absolute value of the Pearson correlation coefficient between the corresponding two factor time series. This graph is drawn using a visualization tool, where the thickness or color intensity of the lines reflects the magnitude of the absolute value of the correlation coefficient, thus visually demonstrating the covariance relationship between environmental factors. The input is the multi-source environmental parameter time series table. For any pair of environmental factors F... i and F j Calculate the effect of controlling for all other environmental factors (i.e., excluding other factors from affecting F). i and F j After the linear effect of the combined influence (parts), F i and F j The partial correlation coefficients between environmental factors are calculated using a method based on the inverse of the covariance matrix (also known as the precision matrix). This involves calculating the ratios of corresponding elements in the precision matrix to obtain the partial correlation coefficients. After calculating the partial correlation coefficients between all environmental factors, an N×N symmetric matrix is ​​constructed, where N is the number of environmental factors. The element in the i-th row and j-th column (and the j-th row and i-th column) of this matrix represents factor F. i and F jThe partial correlation coefficients between environmental factors are analyzed. This matrix reflects the strength of the direct association between environmental factors after removing indirect influences. A threshold is set for the partial correlation coefficients; for example, an absolute value greater than 0.7 indicates a significant direct interaction. If the partial correlation coefficient is greater than the threshold (and positive), a synergistic effect is considered, meaning the factors change in the same direction, and the combined effect is greater than the sum of their individual effects. If the partial correlation coefficient is less than the negative threshold (and negative), an antagonistic effect is considered, meaning the factors change in opposite directions, and each weakens the other's impact on sediment stability. Based on these judgments, factor combinations with significant synergistic or antagonistic effects (which can be pairwise or multi-factor combinations) are identified, and each combination is assigned an interaction pattern type label. These factor combinations with specific interaction patterns and their corresponding pattern types are recorded in a table, forming a factor interaction pattern table. For example, a significant positive biased correlation was identified between rainfall and inflow, which was labeled as a synergistic enhancement model; a significant positive biased correlation was also identified between water level change rate and pore water pressure change rate, which was labeled as a synergistic enhancement model. For each interaction model identified in the factor interaction model table (e.g., the combination of factors A and B under the synergistic enhancement model), a composite time series representing the combined effect was first constructed according to its model type. For example, for synergistic enhancement, the standardized time series of factors A and B can be added together to obtain the composite time series S = standardized(A) + standardized(B). Then, the cross-correlation coefficient between this composite time series S and the resistivity time series of each specific depth layer (e.g., 1 cm, 2 cm...180 cm depth) in the sediment resistivity stratification map was calculated, covering a lag time range of 0 to 72 hours. For each depth layer, the lag time τ that maximizes the absolute value of the cross-correlation coefficient was found. This τ is the optimal time lag for the combined effect of the factor combination on the resistivity of that depth layer. The factor combination under each interaction model, the depth layer it affects, and the corresponding optimal time lag are recorded in a table to form a time lag interaction matrix. For example, the time-delay interaction matrix shows that the synergistic effect of rainfall and inflow is strongest in the 0-30 cm sediment depth layer, with an optimal time delay of 6 hours; the synergistic effect of water level change rate and pore water pressure change rate is strongest in the 50-80 cm sediment depth layer, with an optimal time delay of 12 hours. For each factor combination, depth of influence, and optimal time delay τ recorded in the time-delay interaction matrix, the numerical characteristics of these factor combinations at time τ before the occurrence of the historical instability event are analyzed retrospectively. For example, if the historical landslide event occurred at time T, and the time-delay interaction matrix shows that factor combination (A, B) has a significant impact at depth D with an optimal time delay of τ, then the values ​​of factors A and B at time T-τ are analyzed.Through statistical analysis or rule mining methods, characteristic numerical ranges or relative relationships of these factor combinations are identified at time τ before the occurrence of an unstable event. For synergistic enhancement patterns, it is identified that when factors A and B simultaneously exceed a certain threshold at time T-τ (e.g., A>A0 and B>B0), the instability risk increases significantly. For antagonistic mitigation patterns, it is identified that when one factor changes significantly while the other fails to provide sufficient counterbalancing (e.g., C increases significantly while D does not decrease significantly), the risk increases. The critical numerical conditions or thresholds for each identified factor combination at a specific depth and time lag are recorded in a table, forming a conditional threshold table. For example, the conditional threshold table might specify that when rainfall (6 hours ago) > 20 mm / h and inflow rate change (6 hours ago) > 10 m / s / min, the risk in the shallow layer increases; and when water level change rate (12 hours ago) > 0.1 m / h and pore water pressure (12 hours ago) > 0.05 MPa, the risk in the middle layer increases. The distribution of critical conditions recorded in the condition threshold table at different depths is analyzed. For example, some critical conditions are only effective within a specific depth range. Simultaneously, by combining the influence intensity and optimal time delay information of single factors in the factor-depth response matrix or composite factors in the time-delay interaction matrix at different depths, the significance or sensitivity differences of different interaction modes and their critical conditions in the vertical direction of the sediment are identified. For example, the critical condition of a certain interaction mode may be lower in shallow layers but requires a higher value to trigger in deeper layers; or the influence intensity of a certain interaction mode may peak in a specific depth range. These depth-related interaction characteristics and critical condition information are integrated and visualized to form a depth interaction distribution map. This map can be a table listing key interaction modes, corresponding critical conditions, and influence intensity assessments for different depth ranges (e.g., 0-30cm, 30-80cm, 80-180cm); or a diagram showing the activity or importance of different interaction modes along the depth axis. By integrating all the information from the depth interaction distribution map, the factor interaction mode table, and the condition threshold table, a complete and interpretable model or map, namely the factor interaction effect map, is constructed. This atlas clearly demonstrates: which environmental factors interact (from the factor interaction pattern table); the type of interaction (synergistic or antagonistic); the key critical conditions or threshold combinations that trigger instability signs (from the condition threshold table); and the characteristics and scope of these interactions and critical conditions at different depths of sediment (from the depth interaction distribution map). For example, the factor interaction effect map can be a flowchart or network diagram, where nodes represent environmental factors, edges represent interaction relationships, and the interaction type and critical conditions are labeled on the edges. Different edges or nodes will correspond to different colors or styles to indicate their significance at different depths. This atlas provides a summary and quantitative description of how environmental factors complexly affect sediment stability.

[0054] Preferably, the construction of the sediment layer structure model in step S3 includes:

[0055] Based on the resistivity spectral diagram and the stability influencing factor spectrum of the sediment, the physical parameters of the sediment are inverted to obtain the physical property profile of the sediment.

[0056] Based on the resistivity tomography of the sediment, the resistivity peak value of the sediment physical property profile was detected to obtain the candidate point set of the layer interface.

[0057] The candidate point set of the layer interface is clustered and layered according to its characteristics to obtain a physical characteristic layering scheme;

[0058] Identify the distribution map of the interlayer transition zone of the sediment physical property profile based on the physical property stratification scheme;

[0059] Identify key weak layers in the interlayer transition zone distribution map;

[0060] Underwater terrain fusion processing is performed based on a layering scheme according to physical characteristics to obtain a three-dimensional layered terrain model;

[0061] Layered parameter space interpolation is performed on the three-dimensional layered terrain model to obtain the layered physical parameter field;

[0062] A sediment layer structure model is generated based on the layered physical parameter field and the key weak layer location map.

[0063] In this embodiment of the invention, the inputs are a sediment resistivity strata (containing resistivity data at various spatial locations, depths, and times) and a sediment stability influencing factor spectrum (containing the influence weights of environmental factors on sediment stability and critical conditions). For each spatial location and time point in the sediment resistivity strata, its vertical resistivity profile is extracted. Based on the improved Archie formula, porosity parameters are inverted using the resistivity data. The improved Archie formula uses ρ = a × ρ_w × φ. -m The formula is: ρ = ρ_w / ρ_s, where ρ is the resistivity of the sediment (in ohms·m), ρ_w is the resistivity of pore water (estimated by conductivity or total dissolved solids (TDS) in water quality parameters, in ohms·m), φ is the porosity of the sediment (dimensionless), a is the lithology coefficient (ranging from 0.5 to 2.5, with 1 for pure sediment), and m is the cementation index (ranging from 1.3 to 2.5, with 2 for loose sediment). The values ​​of a and m are determined based on the sediment type (preliminarily determined through historical sedimentary data or resistivity profile morphology). Based on the inverted porosity φ, the water content w and dry density ρ_d are calculated using empirical formulas. Water content w = φ × ρ_w / ρ_s, where ρ_s is the density of solid particles (ranging from 2.65 to 2.75 g / cm³ for sediment). Dry density ρ_d = (1-φ) × ρ_s. Total density ρ_bulk = ρ_d + φ × ρ_w. For cohesion c and internal friction angle φ... iThe estimation is based on the type of sediment and water content. For example, for silty clay, the cohesion c is negatively correlated with the water content w, and the internal friction angle φ... i The strength is negatively correlated with the water content (w), and the specific relationship is represented by a piecewise linear function or exponential function, with parameters calibrated based on historical site test data (such as vane shear tests). Simultaneously, the estimated cohesion and internal friction angle are corrected by referring to the influence weights and critical conditions of environmental factors (such as pore water pressure and water level change rate) on sediment strength in the sediment stability influence factor spectrum, especially when environmental factors approach critical values, reducing the estimated strength parameter values. The density, water content, porosity, cohesion, and internal friction angle values ​​corresponding to each spatial location, time point, and depth are recorded to form a sediment physical property profile dataset. The input is the vertical resistivity profile data from the sediment resistivity stratification diagram. For each spatial location and time point, the first derivative dρ / dz with respect to depth is calculated for the resistivity-depth profile curve at that location. The derivative is calculated using the central difference method, i.e., (ρ(z+Δz)-ρ(z-Δz)) / (2Δz), where Δz is taken as 1 cm. Before calculation, the original resistivity profile was smoothed using a Savitzky-Golay filter with a window length of 5 cm to reduce noise. Local maxima and minima (i.e., zero-crossing points of the second derivative) in the dρ / dz curve were detected. The depths corresponding to these local extrema were considered as the locations with the largest resistivity changes, typically corresponding to layer interfaces where the sediment material properties or water content change significantly. A threshold for the absolute value of the derivative was set, for example, |dρ / dz|>0.05 ohm·m / cm, and only extrema with absolute derivative values ​​exceeding this threshold were retained to exclude minor fluctuations. The depths corresponding to these local extrema exceeding the threshold were recorded as a candidate set of layer interfaces for that spatial location and time point.

[0064] The input consists of a set of candidate points for the interlayer boundary and profile data of the sediment physical properties. For each spatial location and time point, the vertical depth range of that location is divided into several preliminary segments using the candidate points for the interlayer boundary. For each segment, the average or median value of the corresponding sediment physical property profile data (density, water content, porosity, cohesion, internal friction angle) is extracted as the feature vector for that segment. The feature vectors of all spatial locations, all time points, and all segments are aggregated. K-means clustering is used to perform cluster analysis on these feature vectors. The number of clusters K is set to 3, representing the main sediment types such as silt, silty clay, and silt, and this number is set based on historical geological exploration data. The clustering results assign each segment to a specific physical property category (i.e., a cluster). Segments of the same category that are vertically adjacent are merged to form layers with relatively uniform physical properties. For each spatial location and time point, the vertical stratification structure is determined based on the clustering results, including the thickness, burial depth, and physical property category of each layer, forming a physical property stratification scheme. The inputs are the physical property stratification scheme and sediment physical property profile data. For each layer interface (i.e., the boundary between two adjacent layers of different physical property categories) determined in the physical property stratification scheme, a 10 cm depth range is extended above and below this interface as a check window. Within this window, the gradient changes of the sediment physical property profile data (e.g., water content or density) are analyzed. If the absolute value of the gradient of any physical parameter (e.g., water content) within the window exceeds 20% of the average parameter difference between the two main layers above and below the interface (e.g., |dw / dz|>0.2×|w_upper_avg-w_lower_avg|), or the absolute value of the density gradient exceeds 0.03 g / cm³ / cm, then the region is considered an interlayer transition zone. The starting and ending depths of the transition zone are determined, i.e., the depth range that satisfies the gradient threshold condition. The spatial location, time point, depth range of the identified interlayer transition zone, and corresponding physical property changes are recorded to form an interlayer transition zone distribution map. The inputs are the interlayer transition zone distribution map, sediment physical property profile data, and sediment stability influencing factor spectrum. First, the average cohesion *c* and internal friction angle *φ* within the interlayer transition zone are extracted from the sediment physical property profile data. i Based on a preset strength threshold (e.g., cohesion c < 3 kPa or internal friction angle φ), i(<15 degrees) Low-intensity transition zones are marked. Then, referring to the sediment stability influence factor spectrum, depth intervals particularly sensitive to changes in pore water pressure or stress are identified. Low-intensity transition zones are overlaid with depth intervals sensitive to environmental changes. If an interlayer transition zone simultaneously exhibits low intensity and high sensitivity to environmental changes (e.g., absolute correlation with pore water pressure >0.7), it is marked as a critical weak layer. A critical weak layer map is formed by recording the depth range and weakness characteristics (low intensity, high sensitivity) of each spatial location, time point, and identified critical weak layers. Inputs are a physical property stratification scheme and a digital elevation model (DEM) of the reservoir bottom. The DEM of the reservoir bottom is acquired using a high-precision single-beam or multi-beam echo sounder system with a horizontal resolution of 1 m × 1 m and a vertical accuracy of 0.1 m. The physical property stratification scheme provides the interlayer interface depth (relative to the sediment surface) in the vertical direction at each monitoring location (discrete point). It is assumed that the interlayer interfaces of the same physical property category are continuous and smoothly varied in the horizontal direction and are approximately parallel to the bottom topography of the sediment. For each grid point (x, y) on the DEM, first obtain its corresponding bottom elevation Z_bottom(x, y) of the sediment. Then, find the three closest monitoring locations to this grid point and obtain the physical property stratification schemes for these three locations. Based on the relative depth information of the layer interfaces at these three locations and their horizontal distances from the (x, y) point, use the inverse distance weighting method or Kriging interpolation method to estimate the depth d of the layer interface for each physical property category at the (x, y) point relative to the sediment surface. i (x, y). The position of this interface at absolute elevation is Z. i (x,y)=Z_bottom(x,y)-d i(x, y). Interpolation calculations are performed on all layer interfaces and all DEM grid points to obtain the three-dimensional spatial boundary (i.e., layer interface elevation distribution) of each physical property category layer within the entire monitoring area. These boundaries collectively constitute a three-dimensional layered topographic model of the sediment, which describes the layered geometry of the sediment in space. The inputs are the three-dimensional layered topographic model and the sediment physical property profile dataset. The three-dimensional layered topographic model divides the sediment volume into several three-dimensional regions with specific physical property categories. For each physical property category, the physical parameter values ​​(density, water content, porosity, cohesion, internal friction angle) of all layers belonging to that category in the sediment physical property profile dataset are collected. These parameter values ​​are discrete spatial point data (located at a specific depth at the monitoring location). Using the three-dimensional kriging interpolation method, these discrete point data are used to interpolate the parameters of the entire three-dimensional region corresponding to the physical property category in the three-dimensional layered topographic model. The interpolation calculation estimates the corresponding physical parameter value for each element (e.g., finite element mesh element) in the three-dimensional layered topographic model. The interpolation process considers the spatial correlation of the parameters. The resulting layered physical parameter field is a three-dimensional dataset where each spatial location (or mesh cell) is assigned a set of spatially interpolated physical parameter values, including density, water content, porosity, cohesion, and internal friction angle, corresponding to its corresponding physical property category. The inputs are the layered physical parameter field and a critical weak layer map. The layered physical parameter field already provides the three-dimensional physical property distribution of the main sediment mass. Information from the critical weak layer map is overlaid onto this three-dimensional model. The critical weak layer map marks weak regions within specific spatial locations and depth ranges. When generating the sediment layered structure model for subsequent finite element analysis, the volumetric elements corresponding to these critical weak layers are specifically marked. They can be assigned modified, lower strength parameter values ​​(cohesion, internal friction angle) based on their degree of weakness (e.g., assessment of weak features from the critical weak layer map), or these regions can be locally refined during finite element mesh generation. The final generated sediment layer structure model is a digital representation that includes the three-dimensional geometry of the sediment, its internal layer structure, the physical parameters of spatial variations at each location in each layer, and clearly identifies and highlights key weak layers, providing accurate geometric and material property inputs for subsequent stress field calculations.

[0065] Preferably, the calculation of the foundation stress field in step S3 includes:

[0066] Stress-strain relationship calculations were performed on the sediment layer structure model to obtain the set of sediment mechanical response parameters.

[0067] The sediment layer structure model is discretized into a finite element mesh with a horizontal mesh size of 1-5 meters and a vertical mesh size of 5-20 centimeters. The mesh density is adjusted with depth, and the mesh is fined at key layers. Then, each mesh element is assigned corresponding mechanical response parameters according to the sediment mechanical response parameter set, and finally the sediment finite element mesh model is obtained.

[0068] Calculate the initial stress distribution diagram of the finite element mesh model of the bottom mud.

[0069] In this embodiment of the invention, the input is a sediment layer structure model, which provides the geometric layering information of the sediment in space and the physical parameters (density ρ_bulk, porosity φ, water content w, cohesion c, internal friction angle φ) of each layer at each location. i Based on these physical parameters, the stress-strain constitutive relationship of the sediment material is defined. An elastoplastic constitutive model is adopted, such as a model based on the Mohr-Coulomb yield criterion, which follows Hooke's law in the elastic stage and undergoes plastic deformation or failure after reaching the yield surface. The elastic modulus E and Poisson's ratio ν are calculated and determined. The elastic modulus E can be estimated empirically using the sediment density ρ_bulk and porosity φ, for example, E = k × (ρ_bulk). 2 / φ p Where k and p are empirical coefficients, their values ​​determined based on the type of sediment (determined through a stratification scheme using physical properties) and historical site test data. Poisson's ratio ν for saturated soft soil is typically taken as a value close to 0.5, for example, 0.45. Cohesion c and internal friction angle φ... i These parameters (E, ν, c, φ) are directly taken from the corresponding physical parameter values ​​in the sediment layering structure model. i These parameters together constitute a description of the mechanical response characteristics of the sediment at different locations and under different conditions. Linking these parameter sets with their spatial locations in the sediment layering structure model forms the sediment mechanical response parameter set, which is a spatially distributed parameter field.

[0070] The input consists of a sediment layering structure model and a sediment mechanical response parameter set. Using a 3D finite element preprocessing tool, the sediment layering structure model is discretized using a structured or unstructured mesh generation algorithm based on its geometric boundaries. Three-dimensional solid elements (e.g., hexahedral or tetrahedral elements) are generated. In the horizontal direction, the side length of the mesh element is set between 1 and 5 meters, with the specific value determined based on the size of the reservoir area and the required computational accuracy. In the vertical direction, the height of the mesh element is set between 5 and 20 centimeters to precisely capture stress changes in the vertical direction. Based on the key weak layers identified in the sediment layering structure model (from the sub-step of building the sediment layering structure model in S3), the mesh in these areas is locally refined, reducing the mesh size; for example, the vertical mesh size can be reduced to 2 centimeters. After mesh generation, all generated finite element mesh elements are traversed. For each mesh element, its spatial location (e.g., element center coordinates) is determined. Based on the layer and physical property category to which this location belongs in the sediment layering structure model, the corresponding physical parameters (E, ν, c, φ) in the sediment mechanical response parameter set are found. i These parameters are then assigned as material properties to the mesh element. The resulting finite element mesh model of the sediment is a digital representation that includes node coordinates, element connectivity, and the material properties (mechanical response parameters) of each element.

[0071] The input is a finite element mesh model of the sediment. Based on this model, the initial stress state of the sediment under its own weight and water pressure is calculated using a finite element solver. The density ρ_bulk of the sediment material is applied as a volume load to each mesh element to simulate gravity. Water pressure is applied to the sediment surface, and the water pressure value is calculated based on the current reservoir water level, increasing linearly with depth (p = ρ_w × g × h, where p is the water pressure (unit: Pascal), ρ_w is the density of water (unit: kg / m³), and g is the acceleration due to gravity (taken as 9.81 m / s²). 2 (where h is the depth from the water surface (in meters)). Set boundary conditions, such as fixing the vertical displacement at the bottom boundary of the sediment and constraining the horizontal displacement at the lateral boundary (or applying boundary conditions related to horizontal stress, such as the K0 condition). Solve the static equilibrium equation [K]{u}={F}, where [K] is the global stiffness matrix, {u} is the nodal displacement vector, and {F} is the external load vector (including body load and surface load). After calculating the displacements of all nodes, calculate the stress components (normal stress σ) of each element or node based on the strain-displacement relationship of the element and the stress-strain relationship of the material. x ,σ γ ,σ2 and shear stress τ xγ ,τ γ2 ,τ 2xThese stresses are the total stresses. The calculated stress component values ​​and their spatial location information for each node or unit of the sediment are stored and visualized to form an initial stress distribution map of the sediment.

[0072] Preferably, the pore water pressure calculation in step S3 includes:

[0073] Dynamic permeability coefficient correction was performed on the resistivity stratification diagram and the stability influencing factor spectrum of the sediment to obtain the dynamic permeability coefficient field.

[0074] Based on the dynamic permeability coefficient field, a multi-source water pressure superposition calculation is performed to obtain a pressure component analysis table;

[0075] Based on the pressure component analysis table, the interface pressure jump is identified, and the pressure jump location map is obtained.

[0076] The pressure wave velocity distribution diagram was obtained by analyzing the pressure component analysis table.

[0077] Based on the pressure wave velocity distribution map and the pressure jump location map, the critical instability pressure is predicted, and the pore pressure risk map is obtained.

[0078] A pore water pressure distribution map is generated based on the initial stress distribution map and the pore pressure risk map.

[0079] The effective stress distribution diagram is calculated based on the pore water pressure distribution diagram and the initial stress distribution diagram.

[0080] In this embodiment of the invention, the inputs are a sediment resistivity stratification map (containing resistivity data at various spatial locations, depths, and times) and a sediment stability influence factor spectrum (containing the influence weights of environmental factors on sediment stability and critical conditions). First, based on the resistivity values ​​in the sediment resistivity stratification map, the sediment porosity φ is calculated using an empirical relationship based on the Archie formula, for example, φ = (a × ρ_w / ρ)^(1 / m), where ρ is resistivity, ρ_w is pore water resistivity, and a and m are empirical coefficients whose values ​​are set according to the sediment type (determined through the sub-step of constructing the sediment stratification structure model in step S3). Then, based on the calculated porosity φ and the sediment type, the initial permeability coefficient k0 of the sediment is estimated using the Kozeny-Carman equation: k0 = C × φ 3 / (1-φ) 2Where C is a constant related to particle shape and size distribution, the value of which is determined according to the sediment type. This k0 is a permeability coefficient based on static physical properties. Next, referring to the sediment stability influence factor spectrum, environmental factors (e.g., water level change rate, rainfall intensity) that have a significant impact on pore water pressure changes, along with their critical thresholds and influence weights, are identified. When these environmental factors approach or exceed their critical thresholds, the permeability of the sediment changes dynamically (e.g., rapid water level drops lead to surface cracking and increased permeability, while sustained high water pressure causes particle migration and pore blockage, reducing permeability). A dynamic correction factor F_dynamic is calculated based on the deviation of the current value of the environmental factor from its critical threshold and its influence weight. For example, F_dynamic = 1 + W_factor1 × f1(V_factor1) + W_factor2 × f2(V_factor2) + ..., where W is the influence weight, and f is a function of the environmental factor value V relative to the critical threshold (e.g., the greater the exceedance of the threshold, the larger the function value). The dynamic permeability coefficient k = k0 × F_dynamic. The dynamic permeability coefficient k calculated at each spatial location, depth, and time point is stored to form a dynamic permeability coefficient field, which is a four-dimensional (x,y,z,t) parameter field.

[0081] The input consists of a dynamic permeability coefficient field and a time series table of multi-source environmental parameters (including data on water level changes, rainfall, etc.). A three-dimensional transient seepage finite element model is used, discretizing the sediment volume into a mesh compatible with stress field calculations. The dynamic permeability coefficient field is input as a material property into the seepage model. Based on the time series table of multi-source environmental parameters, time-varying boundary conditions are applied to the boundaries of the seepage model: a time-varying head boundary (head h = water level elevation) related to the reservoir water level is applied to the sediment surface (bottom of the reservoir); a flow boundary related to rainfall intensity is applied to infiltration boundaries (such as rainfall areas); and a zero flow boundary is applied to impermeable boundaries (such as bedrock interfaces). The transient seepage equation is then solved. Where ρ_w is the density of water, S_s is the storage ratio (related to the compressibility and porosity of the sediment), h is the hydraulic head, k is the dynamic permeability coefficient, and Q is the source-sink term (such as rainfall infiltration). The hydraulic head h(x,y,z,t) at each point inside the sediment as a function of time is obtained. This hydraulic head is then converted into pore water pressure u(x,y,z,t) = ρ_w × g × h(x,y,z,t), where g is the acceleration due to gravity. The calculated pore water pressure field is analyzed. It can be decomposed into hydrostatic pressure components (pressure caused only by the current water level) and excess hydrostatic pressure components (additional pressure caused by transient effects such as seepage, water level changes, and rainfall). The total pore water pressure value at each spatial location, depth, and time point, along with (optionally) its decomposed components, is recorded in a table, forming a pressure component analysis table.

[0082] The input is a pressure component analysis table (containing pore water pressure u(x,y,z,t) data). For each spatial location and time point, analyze the pore water pressure distribution curve u(z) along the vertical depth direction. Calculate the gradient of pore water pressure with respect to depth. Identification There is a significant abrupt change in the curve (i.e., the second derivative). The depth locations where local extrema or sign variations exist. These locations typically correspond to sediment interfaces or areas where permeability coefficients change significantly. A threshold is set for abrupt changes in the pressure gradient, for example... The rate of change exceeds a certain percentage threshold, or the absolute value of the second derivative exceeds a certain threshold. The depth locations where significant pressure jumps are identified, along with their corresponding spatial coordinates and timestamps, are recorded to form a pressure jump location map. This map marks the interface locations within the sediment where abnormal seepage or stress concentration may occur.

[0083] The input is a pressure component analysis table (containing pore water pressure u(x,y,z,t) data). Identify typical transient events occurring in the multi-source environmental parameter time series, such as a rapid drop in water level or a significant rainfall event. Track the propagation of pore water pressure disturbances caused by these events within the sediment in the pressure component analysis table. For a specific spatial location, analyze the changes in pore water pressure time series at different depths, identifying the arrival times of pressure peaks or wavefronts at different depths. Calculate the average velocity of the pressure wave propagating from one depth to another, i.e., Δz / Δt, where Δz is the depth difference and Δt is the time required for the pressure wave to propagate. Repeat this process to calculate the pressure wave velocity at different spatial locations and depth ranges (e.g., shallow, middle, and deep layers) under different environmental events. Record the calculated pressure wave velocities and their corresponding spatial locations, depth ranges, and times to form a pressure wave velocity distribution map. This map reflects the speed of pore pressure changes propagating within the sediment and is closely related to the sediment's permeability and compressibility.

[0084] The input includes a pressure wave velocity distribution map, a pressure jump location map, and a sediment stability influencing factor spectrum. Analyze the pressure wave velocity distribution map to identify areas of abnormal pressure wave velocity (e.g., excessively high velocity indicates a dominant seepage channel or high permeability area, while excessively low velocity indicates reduced permeability or stress concentration). Analyze the pressure jump location map to identify interfaces with significant and persistent pressure jumps. Combine this with information from the sediment stability influencing factor spectrum regarding the instability caused by pore water pressure reaching a critical value or pore water pressure ratio reaching a critical value (e.g., u / σ_z approaching or exceeding 1). Compare the current pore water pressure u(x,y,z,t) with the total stress σ_z(x,y,z) at that location (obtained from the initial stress distribution map or considering subsequent stress changes) to calculate the pore water pressure ratio r_u = u / σ_z. Identify areas where r_u approaches or exceeds the critical threshold. Perform an overlay analysis of areas with high r_u, interfaces with significant pressure jumps, and areas of abnormal pressure wave velocity. If a region simultaneously meets multiple anomalous conditions (high r_u, pressure jump, wave velocity anomaly), then the pore water pressure state of that region is considered close to a critical instability state, with a high risk level. Based on the r_u value, pressure jump amplitude, and wave velocity anomaly degree, risk scores or classifications are performed on different regions of the sediment. The spatial distribution of the risk scores or levels is represented to form a pore pressure risk map, which visually displays the regions within the sediment with a high risk of instability due to abnormal pore water pressure.

[0085] The inputs are an initial stress distribution map (containing the total stress σ(x,y,z)) and a pore pressure risk map. While the preceding steps have already calculated the pore water pressure field u(x,y,z,t) (in the pressure component analysis table), this step implies that the total stress and risk information will need to be combined to "generate" the final pore water pressure distribution map. This means using the calculated pore water pressure field u(x,y,z) (at the current time point) and validating, correcting, or highlighting pore water pressure values ​​in high-risk areas based on information from the pore pressure risk map. For example, if the calculated pore water pressure shows a high-risk area in the risk map, that area will be specially marked or color-coded in the final generated pore water pressure distribution map. Alternatively, this step aims to ultimately output a clear and easily understandable spatial distribution map u(x,y,z) of sediment pore water pressure at the current moment, supplemented with risk information. Considering that subsequent steps require accurate pore water pressure values ​​to calculate effective stress, the most reasonable explanation is that this step outputs a precise spatial distribution map of pore water pressure at the current moment, u(x,y,z), which has been verified or confirmed by the aforementioned analysis (including risk assessment).

[0086] The inputs are a pore water pressure distribution map (containing the pore water pressure u(x,y,z) at each point in the sediment at the current moment) and an initial stress distribution map (containing the total stress σ(x,y,z) at each point in the sediment). Based on Terzaghi's effective stress principle, the effective stress tensor σ'(x,y,z) is calculated for each point inside the sediment. The normal stress components (σ...) x ,σ γ The effective stress (σ2) is obtained by subtracting the pore water pressure from the total normal stress: σ' x =σ x -u,σ' γ =σ γ -u, σ'2=σ2-u. Shear stress components (τ) xγ ,τ γ2 ,τ 2x Unaffected by pore water pressure, its effective stress equals the total shear stress: τ' xγ =τ xγ ,τ' γ2 =τ γ2 ,τ' 2x =τ 2x In the calculations, all stress components and pore water pressures must use a consistent sign convention (e.g., compressive stress is positive or negative). The calculated effective stress tensor values ​​at each point in the sediment will then be used to determine the effective stress tensor components. The spatial location information of the sediment is stored to form an effective stress distribution map. This map reflects the actual stress state borne by the sediment skeleton and is a key parameter for judging the strength and deformation of the sediment.

[0087] Preferably, the stress gradient anomaly calculation in step S3 includes:

[0088] Calculate the multi-scale gradient tensor set of the effective stress distribution map;

[0089] Calculate the stress curl distribution map based on the multi-scale gradient tensor set;

[0090] Critical gradient ratio map calculated based on multi-scale gradient tensor set and sediment stratification model;

[0091] Gradient time evolution processing is performed on the multi-scale gradient tensor set to obtain a gradient evolution trend map;

[0092] A stress gradient anomaly region map is generated based on the critical gradient ratio map, stress curl distribution map, and gradient evolution trend map.

[0093] In this embodiment of the invention, the input is an effective stress distribution map, which contains six independent components of the effective stress tensor σ'(x,y,z,t) of each node or element in the sediment finite element mesh model. Data varying with space and time. For each time point, the spatial gradient of the effective stress distribution map is calculated. Gradient operator. In a three-dimensional Cartesian coordinate system, it is represented as Apply the gradient operator to each component of the effective stress tensor, for example, to calculate the gradient of the effective normal stress σ'x. These gradient components collectively constitute the information of the stress tensor gradient. The calculation employs numerical differentiation methods, such as central difference or finite difference methods based on finite element meshes. To achieve "multi-scale" computation, Gaussian smoothing filters with different standard deviations σ_s are applied to the effective stress component field before calculating the gradient. For example, three-dimensional Gaussian filters with standard deviations of 0.5 m, 1 m, and 2 m are used to convolve the effective stress component field, and then the gradient is calculated on each smoothed field. The standard deviation σ_s represents the smoothing scale; smaller σ_s retains more detail, while larger σ_s focuses on macroscopic trends. The stress component gradient vectors calculated at different scales (e.g., for σ'x, at scale i) are then compared. These are collected together to form a multi-scale gradient tensor set, which contains stress tensor gradient information at each spatial location, time point, and different scales.

[0094] The input is a multi-scale gradient tensor set. Stress curl describes the local rotational or vortex tendency of the stress field and is closely related to the shear deformation and potential failure of a material. Although stress is a tensor, the curl associated with shear stress can be calculated. For example, consider the shear stress vector. Calculate its curl This partial derivative information can be extracted from multi-scale gradient tensor sets (e.g., yes (The y-component). At each selected scale, the three components of the curl vector are calculated based on data from the multi-scale gradient tensor set. The stress curl vector calculated at each spatial location, time point, and scale is stored.

[0095] Typically, the magnitude of the curl vector is of interest.

[0096] As a key component of the stress curl distribution map, it represents the rotational intensity of the local stress field. This map reflects the non-uniformity and vortex degree of the local shear stress field within the sediment.

[0097] The input consists of a multi-scale gradient tensor set and a sediment layering model (which provides the cohesion c and internal friction angle φ at each location). i The critical gradient ratio is used to measure the severity of the stress gradient relative to the material strength. For each point within the sediment, its local shear strength is calculated. in This is the effective normal stress at that point (obtained from the effective stress distribution map). A representative stress gradient modulus is selected from the multi-scale gradient tensor set; for example, the stress tensor gradient norm at a specific scale (e.g., a scale close to the thickness of the critical weak layer). (The Frobenius norm is defined as the square root of the sum of the squares of the components of the tensor) or the maximum shear stress gradient modulus. Calculate the critical gradient ratio or A high R value indicates a large stress gradient in regions with low shear strength. The critical gradient ratio R calculated for each spatial location and time point is stored to form a critical gradient ratio map. This map visually shows which regions within the sediment exhibit abnormally high stress change rates compared to their local strength; these regions represent potential stress concentrations and failure initiation sites.

[0098] The input is a multi-scale gradient tensor set (containing gradient information at different time points). For each spatial location within the sediment and for each selected scale, the variation of the stress gradient tensor (or its modulus and components) with time is analyzed. The rate of change of the stress gradient with time, i.e., the time derivative, is calculated. or The time derivative is calculated using time series differencing methods, such as dividing the difference in gradient values ​​between two consecutive time points by the time interval. To smooth out time series noise, a moving average with a window length of three time points can be applied to the gradient time series before calculating the time derivative. Regions where the gradient increases rapidly over time are identified. The gradient is positive and exceeds a certain threshold. The calculated gradient time derivative or its indicated evolution trend (e.g., rapid increase, stability, rapid decrease) at each spatial location and scale is stored to form a gradient evolution trend map. This map reveals the dynamic changes in the stress gradient within the sediment over time; a rapidly increasing gradient trend is an important precursor to impending instability.

[0099] The input consists of a critical gradient ratio map, a stress curl distribution map (focusing on the curl modulus), and a gradient evolution trend map. By comprehensively analyzing the information from these three maps, regions of abnormal stress gradients are identified. Threshold conditions for judging anomalies are set. For example, a region is marked as abnormal if it simultaneously meets the following conditions: ① the critical gradient ratio R exceeds the threshold R_threshold (e.g., R>0.8); ② the stress curl modulus... Exceeding the threshold Curl_threshold (e.g., Pascals per meter (Pa); ③ The gradient evolution trend indicates that the gradient is increasing rapidly, that is... Exceeding the threshold Rate_threshold (for example, (Pascals / meter / hour). These thresholds are determined based on a historical database of precursory disaster features (from the concept in step S4.1) or engineering experience. Identify continuous regions in three-dimensional space that meet the abnormal conditions. Classify the abnormal regions (e.g., Level 1 anomaly, Level 2 anomaly) or assign an anomaly index based on the degree of condition satisfaction, region size, and the magnitude of the anomaly value. Visualize the spatial location, shape, anomaly level, or index of the identified abnormal regions to form a stress gradient anomaly region map. This map clearly indicates areas with abnormal internal stress states in the sediment that are developing towards instability, serving as a key basis for early warning.

[0100] Preferably, the potential sliding surface stress evolution analysis in step S3 includes:

[0101] Potential sliding surfaces are identified by analyzing the stress gradient anomaly region map and the intensity ratio distribution map, resulting in a potential sliding surface distribution map.

[0102] Calculate the sliding surface safety factor table for the potential sliding surface distribution map;

[0103] Based on the effective stress distribution map, the sliding surface safety factor table, and the stability influence factor spectrum, the spatiotemporal evolution analysis of the stress field was conducted to obtain the stress field spectrum between the bottom mud layers.

[0104] In this embodiment of the invention, the inputs are a stress gradient anomaly region map (marking the location and intensity of stress gradient anomalies) and a strength ratio distribution map (describing the spatial distribution of the ratio of shear stress to shear strength). First, these two maps are superimposed in three-dimensional space. Mesh cells or regions that simultaneously meet the following conditions are identified: their strength ratio is close to or exceeds a preset threshold (e.g., strength ratio > 0.7) and they are located within a stress gradient anomaly region (e.g., anomaly index > 0.5). These regions indicate locations where material strength is relatively low and stress changes drastically, representing potential starting points for failure initiation. Next, a graph-based connected domain analysis method is employed. Mesh cells meeting the above conditions are considered nodes in the graph. If two nodes are spatially adjacent and their connection direction is approximately consistent with the principal shear stress direction of the region (obtained from the effective stress distribution map) (e.g., the angle is less than 30 degrees), a connection edge is established between them. On the constructed connected graph, continuous connected regions of a certain size (e.g., containing more than 100 mesh cells) are searched; these regions are considered potential slip surfaces. To obtain a more accurate geometry of the sliding surface, surface fitting techniques (e.g., least-squares polynomial surface or spline surface fitting) or minimum energy path search algorithms are employed for nodes within the identified potential sliding surface regions. The minimum energy path search seeks paths connecting starting points (e.g., the strongest stress gradient anomaly) and ending points (e.g., free boundaries) in regions with high stress gradients and high strength ratios. The path with the minimum accumulated "energy" (e.g., weights associated with the strength ratio and stress gradient) is considered the optimal sliding surface. The spatial coordinate set of each identified potential sliding surface (e.g., a list of mesh cells constituting the surface or parameters of the fitted surface) and the associated eigenvalues ​​(e.g., average strength ratio, maximum stress gradient anomaly index) are recorded to form a potential sliding surface distribution map.

[0105] The input is a potential slip surface distribution map (containing geometric information of the identified potential slip surfaces). For each potential slip surface identified in the potential slip surface distribution map, its safety factor is calculated using the limit equilibrium method. The sliding body enclosed by the slip surface is divided into a series of vertical or inclined strips (e.g., strips spaced 1 meter apart horizontally). For each strip, the average effective normal stress σ' at its bottom (i.e., on the slip surface) is obtained from the effective stress distribution map. And the average shear stress τ. The effective cohesion c' and effective internal friction angle φ of the bottom material of the strip are obtained from the physical property profile of the sediment (or the layered physical parameter field). i Calculate the anti-slip force of each strip. Among them l i This is the length of the bottom of the strip on the sliding surface. Calculate the sliding force T of each strip. i Typically, it is the weight of the strip itself, W. iThe component along the sliding surface direction also needs to consider external loads (such as additional stress caused by changes in water pressure). The safety factor F_s of the entire sliding body is calculated using a specific method of limit equilibrium (e.g., the simplified Bishop method or the Spencer method). The simplified Bishop method's safety factor formula is: F_s=[∑(c'×l i +(W i -u i ×l i )×tanφ i ) / m γi ] / ∑(W i ×sinα i ), where u i It is the pore water pressure at the bottom of the strip (obtained from the pore water pressure distribution map), α i It is the inclination angle of the bottom sliding surface of the strip, m γi =cosα i +sinα i ×tanφ i / F_s (requires iterative solution). The Spencer method considers the force and moment balance between blocks, resulting in more accurate calculations. A table is created by recording the unique identifier of each potential sliding surface and the calculated safety factor value.

[0106] The inputs are an effective stress distribution map (containing the effective stress field σ'(x,y,z,t) that varies with time), a sliding surface safety factor table (containing the safety factor F_s(t) of each potential sliding surface that varies with time), and a sediment stability influencing factor spectrum (containing the influence mechanism of environmental factors on sediment stability). First, the time series data of the effective stress distribution map are analyzed. For key areas (e.g., areas with abnormal stress gradients, near potential sliding surfaces), the time series of their effective stress components (such as maximum shear stress τ_max' and effective normal stress σ') are extracted, and their time derivatives are calculated. and Identify the locations and types of rapid stress changes. Next, analyze the time-series data of the slip surface safety factor table. For each potential slip surface, plot its F_s versus time curve and calculate the time derivative of the safety factor, dF_s / dt. Pay particular attention to slip surfaces where F_s is close to 1 and dF_s / dt is negative and has a large absolute value; these surfaces are rapidly evolving towards instability. Combine this with the sediment stability influence factor spectrum to correlate the current state of environmental factors (obtained from the multi-source environmental parameter time-series table) with changes in the stress field and safety factor. For example, if the water level is rapidly dropping and dF_s / dt is negative, while the influence factor spectrum indicates that the drop in water level reduces pore water pressure, thereby increasing effective stress (usually increasing the safety factor), further analysis is needed to determine if there are other antagonistic factors or anomalous responses in specific strata causing a decrease in the safety factor, or if the current calculation model is biased. By utilizing current trends in environmental factors and correlation patterns between environmental factors and stress / safety factor changes in historical data, the effective stress variation trends at key locations and the safety factor variation trends of key potential sliding surfaces are predicted over a future period (e.g., the next 24 hours or 72 hours). All this dynamic information—including real-time effective stress distribution at various depths of the sediment, real-time pore water pressure distribution (obtained from pore water pressure calculation in step S3), the current location and morphology of identified potential sliding surfaces, their real-time safety factor values, the time-varying rates of change of key stress parameters and safety factors, and the predicted future trends—is integrated and visualized to generate a sediment interlayer stress field map. This map is a comprehensive and dynamic representation that fully reflects the internal mechanical state of the sediment and its evolution with environmental changes and time, serving as the core basis for catastrophic criticality prediction.

[0107] Preferably, step S4 includes the following steps:

[0108] Step S41: Obtain historical disaster data of reservoir bottom sediment; extract historical disaster features from the historical disaster data of reservoir bottom sediment to obtain a disaster precursor feature library;

[0109] Step S42: Detect stress anomalies based on the interlayer stress field map of the bottom sediment to obtain a table of stress anomaly regions;

[0110] Step S43: Based on the stress field anomaly area table and the disaster precursor feature database, conduct an instability probability assessment to obtain a regional instability risk map;

[0111] Step S44: Calculate the critical index based on the regional instability risk map and the stress field anomaly area table to obtain the bottom sediment catastrophic critical index;

[0112] Step S45: Generate and release early warning information based on the critical index of sediment disaster, and obtain a reservoir sediment disaster early warning report.

[0113] In this embodiment of the invention, records of all sediment-related disaster events since the reservoir's construction are collected, including but not limited to sediment landslides, mudflows, and liquefaction. These records are derived from accident reports, monitoring data, on-site photos, videos, media reports, and relevant research literature from the reservoir management department. For each disaster event, detailed records are kept of its occurrence time, specific location, affected area, disaster type (e.g., shallow sliding, deep sliding, mudflow), triggering factors (e.g., heavy rain, sudden drop in water level, earthquake), and losses caused. If applicable, monitoring data from a period prior to the disaster are retrospectively retrieved, particularly environmental data such as reservoir water level, rainfall, and inflow, as well as (if available) historical monitoring system records of sediment surface displacement and pore water pressure. Using this data, combined with the stress field calculation and potential sliding surface analysis methods in step S3, the sediment stress field state and potential sliding surface characteristics before historical disasters are simulated or inverted. Typical stress field precursor signals before a disaster are extracted, such as: peak shear stress at critical locations, rapid increase in stress gradient, abnormal increase in pore water pressure leading to a significant decrease in effective stress, potential sliding surface safety factor dropping to a critical value (e.g., between 1.1 and 1.3) and continuing to decrease, and rapid decrease in resistivity at specific depths. The extracted typical stress field precursor feature patterns of different types of sediment disasters (e.g., shallow sliding precursors: high intensity ratios concentrated in the shallow layer, rapid increase in shallow stress gradient; deep sliding precursors: continuous decrease in the safety factor of the deep potential sliding surface; mudflow precursors: rapid increase in pore water pressure leading to effective stress approaching zero) and their corresponding environmental triggering conditions and occurrence time windows are structured and stored to form a disaster precursor feature library. This library serves as a benchmark for subsequent anomaly detection and probability assessment.

[0114] The input is the current stress field map of the sediment layers (including real-time effective stress distribution, pore water pressure distribution, potential slip surface location, safety factor, etc.). Based on the stress field precursor signal types and thresholds defined in the disaster precursor feature library, the current stress field map is analyzed. For example, it checks whether there are regions within the sediment where the ratio of shear stress to shear strength exceeds 0.8; whether there are regions where the stress gradient (e.g., the critical gradient ratio calculated based on the stress gradient anomaly in step S3) exceeds the threshold R_threshold and the gradient is rapidly increasing (judged based on the gradient evolution trend map); whether there are regions where the pore water pressure ratio (u / σ_z) exceeds the threshold 0.9; and whether there are slip surfaces in the potential slip surface safety factor table with a safety factor lower than 1.3 and a time derivative dF_s / dt less than -0.01 / hour. For any region or potential slip surface that meets any of these anomaly conditions, it is marked as a stress field anomaly. Record the spatial location (e.g., center coordinates or slip surface ID) of each anomalous region or potential slip surface, the anomaly type (e.g., intensity ratio anomaly, gradient anomaly, pore water pressure anomaly, safety factor anomaly), the anomaly severity (e.g., the magnitude exceeding the threshold), and the rate of change of the anomalous parameters. Summarize this information into a table to form a stress field anomaly region table. This table lists all anomalous signals currently present within the sediment that are similar to historical precursor patterns of catastrophic events.

[0115] The input consists of a table of anomaly regions in the stress field and a database of precursory disaster features. For each anomaly region or potential slip surface in the table, its current stress field characteristics (anomaly type, degree, rate of change) are compared with precursory patterns of different types of disasters in the database. The similarity between the current anomaly features and each historical precursory pattern is calculated using various similarity metrics, such as Euclidean distance, cosine similarity, or rule-based matching. For example, if a shallow region with a high intensity ratio and a rapidly increasing stress gradient is detected, the similarity to a shallow slip disaster precursory pattern is high. Based on the similarity, the probability that the anomaly region or potential slip surface will eventually develop into an actual disaster is assessed. This assessment can be performed using Bayesian networks, support vector machines, or rule-based inference systems. For example, P(disaster | current anomaly feature) = P(current anomaly feature | precursory disaster pattern) × P(precursory disaster pattern) / P(current anomaly feature), where P(precursory disaster pattern) can be determined based on the frequency of historical disaster events. The size and spatial distribution characteristics of the anomaly region are also considered. If multiple interconnected anomalous regions occur simultaneously, or if these anomalous regions are located on critical weak layers, the overall probability of instability increases. For potential slip surfaces, the safety factor F_s is a core indicator; the instability probability exhibits a non-linear negative correlation with F_s and a positive correlation with the absolute value of dF_s / dt. By analyzing historical data, quantitative relationship curves between F_s, dF_s / dt, and the instability probability are established. The instability probability of each spatial location (or grid cell) is calculated, forming a three-dimensional probability distribution field. This probability distribution field is visualized to create a regional instability risk map, where different colors or grayscale values ​​represent the probability of instability in each region.

[0116] The input consists of a regional instability risk map and a table of anomaly areas in the stress field. The sediment catastrophic criticality index is a comprehensive indicator that quantifies the overall instability risk of the current sediment. The calculation of this index considers the following dimensions: ① Instability probability: Extracting the highest probability value or the average probability of high-risk areas from the regional instability risk map. ② Anomaly severity: Extracting the most severe anomaly type, the highest anomaly severity, or the anomaly index from the table of anomaly areas in the stress field. ③ Spatial range: Calculating the proportion of the total volume or area of ​​high-risk areas (e.g., instability probability > 0.6) to the entire monitoring area. ④ Time urgency: Based on the rate of change of key anomaly parameters (such as safety factor, stress gradient), combined with the time required for these parameters to go from the appearance of an anomaly to the occurrence of a catastrophic event from the historical precursor feature database, estimating the time window required to reach the critical state (e.g., the safety factor drops to 1.05). The criticality index can be calculated using a weighted summation or multi-factor product. For example, the catastrophic criticality index = W1 × P_max + W2 × A_max + W3 × S_ratio + W4 / T_window, where P_max is the highest probability of instability, A_max is the highest degree of anomaly, S_ratio is the proportion of high-risk areas, T_window is the estimated time window before reaching the critical state, and W1, W2, W3, and W4 are weighting coefficients, whose values ​​are determined based on the analysis of historical disaster cases. Based on the numerical range of the criticality index, a graded early warning standard is set, for example: index < 20 is normal; 20 ≤ index < 50 is caution; 50 ≤ index < 80 is warning; 80 ≤ index < 100 is danger; and index ≥ 100 is emergency. The calculated value is the sediment catastrophic criticality index, which reflects the overall risk level and urgency of current sediment instability.

[0117] The input is the sediment disaster critical index. Based on the calculated critical index value, determine the current warning level (normal, attention, warning, danger, emergency). Based on the current stress field anomaly area table and regional instability risk map, determine the location of the area with the highest risk, the potential disaster type (e.g., if the anomaly is mainly concentrated in the shallow layer and highly correlated with rainfall, it will be shallow slip; if the deep layer safety factor continues to decline, it will be deep slip), and the estimated impact range (e.g., estimated by the spatial extent of the high-risk area). Based on the time window estimated in step S44, give the expected time period for the disaster to occur. Integrate this information into a structured reservoir sediment disaster warning report. The report should include at least: the current warning level, the warning issuance time, a spatial description of the anomaly area (e.g., area name, latitude and longitude range), main anomaly characteristics (e.g., key slip surface safety factor and rate of change, area with the highest intensity ratio), potential disaster type, estimated impact range, expected occurrence time window, and recommended emergency response measures for this warning level and disaster type (e.g., strengthen patrols, limit the reservoir water level drawdown rate, activate the emergency plan, notify downstream areas, etc.). For warnings, hazards, and emergencies, the system automatically sends alerts to designated reservoir management officials and on-duty personnel via SMS, email, and telephone. Simultaneously, it generates intuitive visual charts and graphs, including marking risk areas on the reservoir map, displaying time-series graphs of critical indices to illustrate risk trends, and showing curves illustrating the safety factor changes of key potential slip surfaces. These visualizations are released along with the alert report, enabling managers to quickly understand the risk situation and take action. The final output is a reservoir sediment disaster early warning report.

[0118] Of particular importance is the specific assessment of the probability of instability:

[0119] Based on the table of anomaly regions in the stress field, identify the regional instability mode table;

[0120] Based on the regional instability pattern table, the historical similarity of the stress field anomaly region table and the disaster precursor feature database is calculated to obtain a historical similar case matching table.

[0121] Based on the stress field anomaly region table and the regional instability mode table, a shear-slip instability safety assessment is conducted to obtain a sliding surface safety factor assessment table.

[0122] Based on the table of anomaly regions in the stress field and the table of regional instability modes, a liquefaction-type instability risk assessment is conducted to obtain a liquefaction risk assessment table.

[0123] Based on the historical similar case matching table, the sliding surface safety factor assessment table and the liquefaction risk assessment table, the single-point instability probability is calculated to obtain the point instability probability table.

[0124] Based on the point instability probability table, spatial correlation analysis was performed on the stress field spectrum of the sediment layer to obtain a risk spatial correlation map.

[0125] A cascading effect assessment was conducted on the point instability probability table and the risk space correlation diagram to obtain a disaster spread risk map.

[0126] By performing time window predictions on the disaster spread risk map, a regional instability risk map is obtained.

[0127] In this embodiment of the invention, based on the combination of abnormal features in the stress field anomaly region table, the potential instability modes (such as shear slip and liquefaction) of each anomaly region are identified using preset rules or classification models, resulting in a regional instability mode table; the features of the current anomaly region are compared with the features of historical cases of the same mode in the disaster precursor feature database using similarity calculation (cosine similarity or weighted Euclidean distance), resulting in a historical similar case matching table; for regions identified as shear slip modes, the safety factor F_s and the rate of change dF_s / dt of their potential slip surface are extracted, resulting in a slip surface safety factor assessment table; for regions identified as liquefaction modes, the pore water pressure ratio r_u and material sensitivity are extracted to assess the liquefaction risk, resulting in a liquefaction risk assessment table; and the results are combined with historical similarity and slip surface safety factor assessment. Based on the numerical assessment and liquefaction risk assessment results, the local probability of instability in each anomalous area (location) is calculated to obtain a location instability probability table. Spatial interpolation (such as co-kriging) and spatial correlation analysis are performed on the location instability probabilities to generate a risk spatial correlation map of the regional instability probability distribution. High-risk areas are identified as potential triggering areas, and their instability is simulated to have cascading effects (such as impact and pore pressure diffusion) on surrounding areas (downstream and adjacent areas). Additional risks are superimposed to obtain a catastrophic diffusion risk map considering diffusion effects. Finally, combining the catastrophic diffusion risk map and the time change rate of key precursor parameters (such as safety factor and stress gradient), historical data is used to predict the time window for reaching the critical instability state, ultimately forming a regional instability risk map that includes spatial distribution and time prediction.

[0128] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0129] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for predicting and analyzing the operational status of a reservoir based on environmental data, characterized in that, Includes the following steps: Step S1: Using a vertical resistivity probe array, low-frequency AC signals are used to collect the original resistivity response sequences of the reservoir bottom sediment at different depths; polarization effect compensation processing is performed on the original resistivity response sequences to obtain the bottom sediment resistivity spectral map. Step S2: Collect reservoir environmental data; perform time-series coupling analysis on the reservoir environmental data and sediment resistivity spectral diagram to obtain the sediment stability influencing factor spectrum; Step S3: Invert the physical parameters of the sediment based on the sediment resistivity tomography and the sediment stability influencing factor spectrum to obtain a sediment physical property profile; detect the resistivity peak of the sediment physical property profile based on the sediment resistivity tomography to obtain a set of candidate points for the layer interface; perform characteristic clustering and stratification on the candidate points for the layer interface to obtain a physical property stratification scheme; identify the distribution map of the interlayer transition zone of the sediment physical property profile based on the physical property stratification scheme. Identify key weak layers in the interlayer transition zone distribution map; Underwater terrain fusion processing is performed based on a physical property layering scheme to obtain a three-dimensional layered terrain model; layered parameter space interpolation is performed on the three-dimensional layered terrain model to obtain a layered physical parameter field; A sediment layer structure model is generated based on the layered physical parameter field and the key weak layer map. The initial stress distribution map was obtained by performing basic stress field calculations on the sediment layer structure model. The pore water pressure is calculated based on the initial stress distribution diagram to obtain the effective stress distribution diagram; Calculate the strength ratio distribution of the bottom sediment based on the effective stress distribution diagram; Stress gradient anomaly is calculated based on the effective stress distribution map to obtain the stress gradient anomaly region map; By analyzing the stress gradient anomaly region map and the intensity ratio distribution map, the stress evolution of the potential sliding surface is obtained, and the stress field spectrum between the bottom mud layers is obtained. Step S4: Based on the stress field spectrum between the bottom sediment layers, assess the probability of reservoir instability and obtain a regional instability risk map; Based on the regional instability risk map, a disaster threshold early warning is issued, resulting in a reservoir bottom sediment disaster early warning report.

2. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Obtain and deploy vertical resistivity probe arrays in key monitoring areas of the reservoir based on the coordinates of the reservoir bottom sediment area and historical sedimentation data. Each probe array contains 10-15 measuring electrodes at different depths with a depth interval of 10-20 cm, and obtain a probe position depth table. Step S12: Send a low-frequency AC signal of 0.1-10Hz to the probe array according to the probe position depth table. The signal strength is 5-20mA. Then perform 5 repeated measurements at each depth point with a sampling rate of 100Hz and each measurement lasts for 30 seconds. Finally, obtain the original resistance response sequence. Step S13: Perform polarization effect compensation on the original resistance response sequence to obtain a set of corrected resistivity values; Step S14: Perform depth interpolation on the set of corrected resistivity values ​​to obtain the depth distribution curve; Step S15: Construct a resistivity stratum map of the bottom sediment based on the depth distribution curve.

3. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, Step S2 includes the following steps: Step S21: Perform environmental data preprocessing on the reservoir environmental data to obtain a time series table of multi-source environmental parameters; Step S22: Perform hysteresis correlation analysis on the sediment resistivity spectral diagram and the time series table of multi-source environmental parameters to obtain the factor-depth response matrix; Step S23: Perform environmental factor interaction analysis on the factor-depth response matrix to obtain the factor interaction effect diagram; Step S24: Extract the critical conditions for stability from the factor interaction effect diagram and the factor-depth response matrix to obtain the spectrum of factors affecting sediment stability.

4. The reservoir operation status prediction and analysis method based on environmental data according to claim 3, characterized in that, Step S2, the analysis of the interaction effects of environmental factors, includes: Construct a factor correlation network graph based on the factor-depth response matrix; Partial correlation analysis was performed on the factor correlation network diagram to obtain the direct influence matrix of the factors; The interaction pattern classification and identification of the direct influence matrix of factors is performed to obtain the factor interaction pattern table; Based on the factor interaction pattern table, the time-delay interaction effect is analyzed to obtain the time-delay interaction matrix; Determine the conditional threshold table based on the time-delay interaction matrix; A deep difference analysis of the factor-depth response matrix was performed based on the conditional threshold table to obtain a deep interaction distribution map; A factor interaction effect diagram is drawn based on the deep interaction distribution diagram, the factor interaction pattern table, and the conditional threshold table.

5. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, Step S3, the calculation of the foundation stress field, includes: Stress-strain relationship calculations were performed on the layered sediment structure model to obtain the set of sediment mechanical response parameters. The sediment layer structure model is discretized into a finite element mesh with a horizontal mesh size of 1-5 meters and a vertical mesh size of 5-20 centimeters. The mesh density is adjusted with depth, and the mesh is fined at key layers. Then, each mesh element is assigned corresponding mechanical response parameters according to the sediment mechanical response parameter set, and finally the sediment finite element mesh model is obtained. Calculate the initial stress distribution diagram of the finite element mesh model of the bottom sediment.

6. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, The pore water pressure calculation in step S3 includes: Dynamic permeability coefficient correction was performed on the resistivity stratification diagram and the stability influencing factor spectrum of the sediment to obtain the dynamic permeability coefficient field. Based on the dynamic permeability coefficient field, a multi-source water pressure superposition calculation is performed to obtain a pressure component analysis table; Based on the pressure component analysis table, the interface pressure jump is identified, and the pressure jump location map is obtained. The pressure wave velocity distribution diagram was obtained by analyzing the pressure component analysis table. Based on the pressure wave velocity distribution map and the pressure jump location map, the critical instability pressure is predicted, and the pore pressure risk map is obtained. A pore water pressure distribution map is generated based on the initial stress distribution map and the pore pressure risk map. The effective stress distribution diagram is calculated based on the pore water pressure distribution diagram and the initial stress distribution diagram.

7. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, The stress gradient anomaly calculation in step S3 includes: Calculate the multi-scale gradient tensor set of the effective stress distribution map; Calculate the stress curl distribution map based on the multi-scale gradient tensor set; Critical gradient ratio map calculated based on multi-scale gradient tensor set and sediment stratification model; Gradient time evolution processing is performed on the multi-scale gradient tensor set to obtain a gradient evolution trend map; A stress gradient anomaly region map is generated based on the critical gradient ratio map, stress curl distribution map, and gradient evolution trend map.

8. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, The analysis of the evolution of potential sliding surface stress in step S3 includes: Potential sliding surfaces are identified by analyzing the stress gradient anomaly region map and the intensity ratio distribution map, resulting in a potential sliding surface distribution map. Calculate the sliding surface safety factor table for the potential sliding surface distribution map; Based on the effective stress distribution map, the sliding surface safety factor table, and the stability influence factor spectrum, the spatiotemporal evolution analysis of the stress field was conducted to obtain the stress field spectrum between the bottom mud layers.

9. The method for predicting and analyzing the operational status of a reservoir based on environmental data according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Obtain historical disaster data of reservoir bottom sediment; extract historical disaster features from the historical disaster data of reservoir bottom sediment to obtain a disaster precursor feature library; Step S42: Detect stress anomalies based on the interlayer stress field map of the bottom sediment to obtain a table of stress anomaly regions; Step S43: Based on the stress field anomaly area table and the disaster precursor feature database, conduct an instability probability assessment to obtain a regional instability risk map; Step S44: Calculate the critical index based on the regional instability risk map and the stress field anomaly area table to obtain the bottom sediment catastrophic critical index; Step S45: Generate and release early warning information based on the critical index of sediment disaster, and obtain a reservoir sediment disaster early warning report.

Images