Reservoir operation state prediction analysis method based on environmental data

Abstract

The invention relates to the technical field of reservoir safety monitoring and early warning, in particular to a reservoir operation state prediction analysis method based on environmental data. The method comprises the following steps: collecting original resistance response sequences of reservoir bottom mud at different depths by adopting low-frequency alternating-current signals through a laid vertical resistivity probe array; performing polarization effect compensation processing on the original resistance response sequence to obtain a sediment resistivity layer spectrogram; collecting reservoir environment data; performing time sequence coupling analysis on the reservoir environment data and the sediment resistivity layer spectrogram to obtain a sediment stability influence factor spectrum; constructing a bottom mud layered structure model according to the bottom mud resistivity layer spectrogram and the bottom mud stability influence factor spectrum; and carrying out basic stress field calculation on the bottom mud layered structure model to obtain an initial stress distribution diagram. Through a full-chain technical route from sediment internal monitoring to multi-factor comprehensive analysis, the early warning timeliness, accuracy and reliability of reservoir sediment disasters are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field of reservoir safety monitoring and early warning, and in particular to a reservoir operation status prediction and analysis method based on environmental data. Background Art

[0002] Traditional monitoring systems generally have the problem of early warning lag. Most monitoring methods can only capture the surface displacement or deformation of the sediment. There is a time difference between the change in the internal stress state of the sediment and the appearance of obvious displacement on the surface, resulting in a short early warning time window, which makes it difficult to provide sufficient preparation time for emergency response. Traditional technologies have limited monitoring capabilities for the internal structure of the sediment, and mainly rely on indirect methods such as surface displacement monitoring and tilt monitoring. There are blind spots in monitoring the physical state and structural changes of deep sediments, and it is impossible to grasp the early signs of internal instability of the sediment in real time. Existing methods usually analyze environmental factors in isolation, ignoring the interaction between multiple environmental factors and their combined impact on sediment stability, resulting in insufficient accuracy of early warning judgments under complex environmental conditions, and it is difficult to establish an accurate correlation model between environmental changes and sediment stability.

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

[0004] Based on this, it is necessary to provide a reservoir operation status prediction and analysis method based on environmental data to solve at least one of the above technical problems.

[0005] To achieve the above objectives, a reservoir operation status prediction and analysis method based on environmental data includes the following steps:

[0006] Step S1: Using a vertical resistivity probe array, a low-frequency AC signal is used to collect the original resistance response sequence of the reservoir sediment at different depths; the original resistance response sequence is subjected to polarization effect compensation processing to obtain a sediment resistivity layer spectrum;

[0007] Step S2: collecting reservoir environmental data; performing time-series coupling analysis on the reservoir environmental data and the sediment resistivity layer spectrum to obtain a sediment stability influencing factor spectrum;

[0008] Step S3: Constructing a sediment layered structure model based on the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum; performing basic stress field calculation on the sediment layered structure model to obtain an initial stress distribution map; performing pore water pressure calculation based on the initial stress distribution map to obtain an effective stress distribution map; calculating a sediment strength ratio distribution map based on the effective stress distribution map; performing stress gradient anomaly calculation based on the effective stress distribution map to obtain a stress gradient anomaly area map; performing potential sliding surface stress evolution analysis on the stress gradient anomaly area map and the strength ratio distribution map to obtain a sediment interlayer stress field map;

[0009] Step S4: Evaluate the probability of reservoir instability based on the inter-layer stress field map of the sediment to obtain a regional instability risk map; perform a critical disaster warning based on the regional instability risk map to obtain a reservoir sediment disaster warning report.

[0010] This method improves raw data quality through high-precision placement and low-frequency measurement. Accurate polarization and temperature compensation remove interfering factors. High-resolution depth and spatial dynamic interpolation constructs resistivity layer spectra reflecting the internal structure of the sediment, enabling precise, dynamic, and nondestructive monitoring of the sediment's physical state. Rigorous preprocessing of multi-source environmental data ensures data quality. Lagged correlation analysis reveals the transmission lag of environmental influences. Partial correlation and interaction pattern analysis delve deeper into the complex interactions between environmental factors. By extracting critical conditions based on historical events, a more precise quantitative correlation model between environmental change and sediment stability is established, enhancing the scientific nature of the predictions. Based on resistivity inversion of physical parameters and identification of weak zones, a three-dimensional model accurately reflects the heterogeneous structure of the sediment. A refined finite element mesh and transient seepage model accurately calculate the internal stress field and pore water pressure, particularly excess pore water pressure. Multidimensional analysis of stress gradients and temporal tracking of their evolution enable the identification of microscopic stress anomalies before macroscopic displacements. Combined with dynamic assessment of the safety factor of potential sliding surfaces, this method provides a comprehensive, dynamic monitoring map of the internal mechanical state of the sediment, significantly accelerating the detection of signs of instability. The establishment of a historical disaster precursor database provides a benchmark for current anomaly detection, and improves the accuracy of instability probability judgment based on multiple model evaluations and historical similarities; the disaster diffusion risk is evaluated by considering spatial correlation and cascade effects, and the time window for risk outbreak is predicted; finally, by calculating the graded critical index and generating a structured early warning report, the complex analysis results are converted into intuitive and timely early warning information, providing an early, comprehensive and scientific decision-making basis for reservoir safety management, and effectively improving the safety of reservoir operation. Therefore, the present invention provides a reservoir operation status prediction and analysis method based on environmental data, which realizes real-time monitoring of sediment microstructure changes by constructing a sediment internal status monitoring architecture based on resistivity layer spectrum; reconstructs the internal stress field distribution of the sediment using the finite element method to capture early signals of stress anomalies; establishes a time-series coupling model of multi-source environmental factors and sediment stability, and comprehensively analyzes the impact of environmental changes on sediment stability. This full-chain technical route from internal sediment monitoring to multi-factor comprehensive analysis has significantly improved the timeliness, accuracy and reliability of early warning of reservoir sediment disasters. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 The figure is a flowchart of the steps of a reservoir operation status prediction and analysis method based on environmental data.

[0012] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0013] The following is a clear and complete description of the technical method of the present invention in conjunction with the accompanying drawings. It is obvious that the embodiments described are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative efforts are within the scope of protection of the present invention.

[0014] In addition, the accompanying drawings are merely schematic illustrations of the present invention and are not necessarily drawn to scale. Identical reference numerals in the figures denote identical or similar parts, and thus repetitive descriptions thereof will be omitted. Some of the block diagrams shown in the accompanying drawings are functional entities that do not necessarily correspond to physically or logically separate entities. These functional entities may be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor and / or microcontroller approaches.

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

[0016] In the embodiment of the present invention, reference Figure 1 FIG. 1 is a flow chart showing the steps of a method for predicting and analyzing the operation status of a reservoir based on environmental data according to the present invention. In this example, the method for predicting and analyzing the operation status of a reservoir based on environmental data includes the following steps:

[0017] Step S1: Using a vertical resistivity probe array, a low-frequency AC signal is used to collect the original resistance response sequence of the reservoir sediment at different depths; the original resistance response sequence is subjected to polarization effect compensation processing to obtain a sediment resistivity layer spectrum;

[0018] In an embodiment of the present invention, a key monitoring area of ​​a reservoir is determined by using high-precision positioning equipment and historical sedimentation data, and an underwater robot is used to deploy a vertical resistivity probe array containing 10 to 15 measuring electrodes (with a depth interval of 10 to 20 cm). The probes are calibrated using a standard resistance block to obtain a probe position depth table containing spatial coordinates and electrode depths. According to the probe position depth table, a multi-channel measuring instrument is used to send a low-frequency AC signal of 0.1 to 10 Hz and 5 to 20 mA to the probe array, and a four-electrode method is used to repeatedly measure each depth point (5 times, 100 Hz sampling, for 30 seconds) to obtain an original resistance response sequence containing a voltage / current time series. The original resistance response sequence is subjected to time domain differential or complex resistivity processing to eliminate the electrode polarization effect, and temperature compensation based on an empirical formula is performed in combination with temperature sensor data. A set of corrected resistivity values ​​was obtained; a vertical depth distribution curve with a depth resolution of 1 cm was generated for the corrected resistivity value set using cubic spline interpolation; finally, three-dimensional kriging spatial interpolation (horizontal resolution 5 m × 5 m) was performed on the vertical curves at different probe positions, and time series data (1 hour resolution) was integrated to construct a sediment resistivity layer spectrum containing three-dimensional information of time, spatial position and depth.

[0019] Step S2: collecting reservoir environmental data; performing time-series coupling analysis on the reservoir environmental data and the sediment resistivity layer spectrum to obtain a sediment stability influencing factor spectrum;

[0020] In the embodiment of the present invention, multi-source environmental data such as reservoir water quality (pH value, dissolved oxygen, turbidity, water temperature), meteorology (precipitation, air temperature, air pressure) and hydrology (water level, inflow, outflow) are collected, time alignment (hourly interval), missing values ​​are filled with the neighboring mean, outliers are removed by median filtering, and Z-score standardization is performed to obtain a multi-source environmental parameter time series table; the lagged correlation coefficient (0 to 72 hours lag) between the time series of each environmental parameter and the resistivity time series of each depth layer of the bottom mud (1 cm resolution) is calculated, and the correlation stability is evaluated through a 7-day sliding window. The residual correlation is analyzed after STL decomposition of the seasonal parameters to obtain a factor-depth response matrix containing the degree of influence, lag time and stability score; based on A factor-depth response matrix is ​​used to construct a factor correlation network, and partial correlation analysis is used 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. The critical conditions or thresholds of key interaction patterns are determined in combination with historical instability events, and the performance differences of these interaction patterns at different depths are analyzed. Finally, a factor interaction effect diagram is drawn to intuitively display the factor interactions, critical conditions, and depth distribution. Finally, the factor interaction effect diagram, factor-depth response matrix, and historical instability event data are integrated to extract the influence weights of each environmental factor (including single factors and combined factors) on sediment stability and the critical threshold interval that triggers the unstable state, forming a factor spectrum affecting sediment stability.

[0021] Step S3: Constructing a sediment layered structure model based on the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum; performing basic stress field calculation on the sediment layered structure model to obtain an initial stress distribution map; performing pore water pressure calculation based on the initial stress distribution map to obtain an effective stress distribution map; calculating a sediment strength ratio distribution map based on the effective stress distribution map; performing stress gradient anomaly calculation based on the effective stress distribution map to obtain a stress gradient anomaly area map; performing potential sliding surface stress evolution analysis on the stress gradient anomaly area map and the strength ratio distribution map to obtain a sediment interlayer stress field map;

[0022] In the embodiment of the present invention, according to the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum, the improved Archie formula (ρ = a × ρ_w × φ -m ) Inverse the sediment porosity, water content and density, and estimate the cohesion c and internal friction angle φ based on the sediment type and water content i, refer to the influence factor spectrum to correct the intensity parameters, and obtain the sediment physical property profile containing the vertical distribution of various physical parameters; the resistivity profile is differentiated to detect the layer peak to identify the layer interface candidate point, and the physical property feature vectors are extracted from the layer segments divided by the candidate points for K-means clustering (K=3) stratification to identify the interlayer transition zone and the key weak layer with low intensity and high sensitivity; the underwater terrain DEM and the physical property stratification scheme are integrated to construct a three-dimensional geometric model of the interface of each layer through spatial interpolation, and the discrete physical parameters are interpolated by three-dimensional Kriging space to obtain the layered physical parameter field, and the geometric, parameter and weak layer information are integrated to generate the sediment layered structure model; the elastic-plastic constitutive model parameters (E, ν, c, φ) are determined according to the layered physical parameters. i ), discretize the layered model into a finite element grid of 1-5 meters horizontally and 5-20 centimeters vertically (weak layers are encrypted to 2 centimeters), assign material properties, and calculate the initial total stress distribution map under the action of self-weight and water pressure; use the dynamic permeability coefficient field (corrected based on resistivity and environmental factors) and the transient seepage model to calculate the pore water pressure field, and deduct the pore water pressure from the total stress to obtain the effective stress distribution map; calculate the multi-scale spatial gradient, shear stress curl, and stress gradient to local strength ratio of the effective stress field, analyze its time evolution trend, and identify the stress gradient anomaly area map; combine the stress gradient anomaly area and the shear stress to shear strength ratio map to identify the location and morphology of the potential sliding surface, and use the limit equilibrium method (such as the simplified Bishop method) to calculate the safety factor of each potential sliding surface; finally, comprehensively analyze the real-time status, change rate and correlation of the effective stress, pore water pressure, potential sliding surface and safety factor with environmental factors, predict future development trends, and generate a stress field map between sediment layers containing all dynamic mechanical state information.

[0023] Step S4: Evaluate the probability of reservoir instability based on the inter-layer stress field map of the sediment to obtain a regional instability risk map; perform a critical disaster warning based on the regional instability risk map to obtain a reservoir sediment disaster warning report.

[0024] In an embodiment of the present invention, data on historical reservoir sediment disasters (landslides, liquefaction, etc.) are collected, the pre-disaster environment and monitoring data are retrospectively analyzed, the pre-disaster stress field characteristics are simulated or inverted, typical stress field precursor patterns (such as high strength ratio, high gradient, high pore pressure ratio, low Fs and decreasing) are extracted, and a disaster precursor feature library is constructed; according to the current inter-layer stress field map of the sediment, the threshold and pattern of the disaster precursor feature library are compared, the stress field abnormal area is detected, the abnormal type, degree and change rate are recorded, and a stress field abnormal area table is obtained; based on the combination of abnormal features, potential instability modes (shear sliding, liquefaction, etc.) are identified, the potential sliding surface safety factor and its change rate are evaluated for the shear sliding mode, the pore water pressure ratio and material sensitivity are evaluated for the liquefaction mode, and the similarity between the current abnormality and the historical precursor case is calculated; the model evaluation results and the historical similarity are integrated to calculate the local instability probability of each abnormal area using a probability model. The regional instability probability distribution is obtained through spatial interpolation; the possible propagation and impact of disasters are simulated and the risk of disaster spread is assessed by considering the mutual influence of stress / pore pressure and potential cascading effects between regions; the spatial distribution of risks and the change rate of key parameters (such as the safety factor change rate) are integrated to predict the time window when critical states are reached, forming a regional instability risk map that includes spatial risk distribution and time window; based on the highest probability of the regional instability risk map, the scope of the high-risk area, the highest degree of anomaly in the stress field anomaly area table, and the predicted time window, a comprehensive sediment disaster critical index is calculated and the warning level (normal, caution, warning, danger, emergency) is delineated; finally, based on the critical index and regional instability risk map information, a reservoir sediment disaster warning report is generated and released through multiple channels, including the warning level, anomaly location, potential disaster type, impact range, expected occurrence time window, and recommended response measures.

[0025] It is particularly important that the polarization effect compensation process is specifically:

[0026] Extracting the phase difference characteristic matrix of the original resistance response sequence;

[0027] Constructing a frequency response curve for the original resistance response sequence to obtain a frequency response parameter set;

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

[0029] Perform temperature effect correction on the frequency response parameter set according to the temperature distribution profile to obtain a temperature corrected resistivity set;

[0030] Eliminate the electrode contact impedance of the temperature-corrected resistivity set according to the phase difference characteristic matrix to obtain a corrected resistivity value set;

[0031] In one embodiment of the present 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 in the sequence (a specific probe array and depth), a fast Fourier transform (FFT) is applied to extract the complex spectrum values ​​of the voltage and current at an 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 the repeated measurements is stored as a phase difference characteristic matrix, which reflects the imaginary impedance characteristics caused by electrode and sediment polarization. Simultaneously, based on the voltage and current spectrum amplitudes extracted by the FFT, the apparent impedance amplitude |Z| = |V(f0)| / |I(f0)| corresponding to the excitation frequency is calculated, and the average amplitude of the repeated measurements is stored as a frequency response parameter set. Next, real-time temperature data at different depths around the measurement point is collected using a temperature sensor integrated in 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, according to the temperature T at the time of measurement, the empirical temperature correction formula for saturated sediment (e.g., ρ 25 =ρ_a / [1+0.02×(T-25)], where ρ_a is the apparent resistivity calculated from the apparent impedance amplitude and the geometric factor. The apparent resistivity in the frequency response parameter set is temperature-corrected to obtain a temperature-corrected resistivity set, which represents the apparent resistivity at standard temperature. Finally, to eliminate the effects of electrode polarization and sediment-induced polarization on the resistivity amplitude, the complex resistivity theory is used to calculate the temperature-corrected apparent resistivity ρ in the temperature-corrected resistivity set based on the phase difference φ in the phase difference characteristic matrix. 25 , calculate the true resistivity value ρ_corrected = ρ 25 ×cos(φ). This operation extracts the real part of the complex impedance and obtains a set of corrected resistivity values ​​that is closer to the resistive properties of the sediment bulk.

[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 reservoir sediment area coordinates 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-10 Hz with a signal strength of 5-20 mA to the probe array according to the probe position depth table. Repeat the measurement five times for each depth point with a sampling rate of 100 Hz and each measurement lasting 30 seconds. Finally, obtain the original resistance response sequence.

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

[0036] Step S14: performing depth interpolation processing on the corrected resistivity value set to obtain a depth distribution curve;

[0037] Step S15: constructing a sediment resistivity layer spectrum according to the depth distribution curve.

[0038] In an embodiment of the present invention, a high-precision differential GPS receiver is first used to accurately locate the reservoir and obtain WGS84 coordinate data for the reservoir sediment area. Simultaneously, geological survey reports from the reservoir construction period, sediment deposition measurement reports from previous years, and historical landslide accident investigation reports are consulted to extract the average thickness, stratification characteristics, distribution of major soil types, and the location of historically unstable areas of the sediment. A comprehensive analysis of the above coordinate data and historical sedimentation data identifies key monitoring areas within the reservoir sediment area, areas 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 sediment at a specified depth. Each probe array is made of a rigid insulating material (such as a high-strength PVC pipe) onto which 12 ring-shaped stainless steel measuring electrodes are fixed. The electrodes are evenly spaced in the vertical direction, and the electrode spacing is strictly set to 15 cm, thereby achieving stratified monitoring of the sediment from below the surface to a depth of 1.8 meters. After the probe array is deployed, the electrode sequence of each probe array is calibrated using a standard resistor block (for example, a standard tank filled with a KCl solution with a known conductivity of 1413 μS / cm). The probe array is immersed in the standard tank in turn, and the resistance value of each electrode pair under standard conditions is measured using a precision resistivity meter (for example, using the four-electrode method). The deviation between the measured value and the theoretical resistance value of the standard resistor block is recorded, and a calibration coefficient table is generated for subsequent data correction. Finally, the unique identifier of each deployed probe array, its precise latitude and longitude coordinates, and the exact depth of each measuring electrode on the array relative to the bottom mud surface are recorded to form a probe position depth table. A dedicated multi-channel resistivity meter is used to connect to the deployed vertical resistivity probe array via an underwater cable. The meter automatically controls the electrode switching based on the electrode depth information recorded in the probe position depth table, and uses the four-electrode method for measurement. For example, using the Wenner arrangement, current is sent through the outer pair of electrodes (A and B), and the inner pair of electrodes (M and N) measures voltage. The measuring instrument outputs a sinusoidal AC constant current signal with a frequency of 5 Hz and an intensity of 10 mA. For different electrode combinations on the probe array (for example, electrodes 1 and 4 are selected as current-transmitting electrodes A and B, and electrodes 2 and 3 are selected as voltage-measuring electrodes M and N to measure the resistivity of the corresponding depth layer; then the electrode combination is switched 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 collects the time series data of the transmitted current and measured voltage at a sampling rate of 100 Hz. The voltage time series collected in each measurement cycle is divided by the current time series to obtain the time series of the original resistance value. The original resistance value time series data 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.Receive the raw resistance response sequence data. Process the voltage and current time series data for each measurement point and each measurement cycle, and analyze the phase difference between the voltage signal and the current signal. Due to the electrode polarization effect, a frequency-dependent impedance component will be generated in the actual measurement. Use the time domain difference method to calculate the instantaneous rate of change of the voltage time series relative to the current time series within 30 seconds of each measurement cycle, and analyze its attenuation characteristics. Alternatively, use the principle of complex resistivity to identify and separate the imaginary resistance component caused by electrode polarization, and extract the true real resistance component, that is, the resistance value after deducting the influence of the polarization effect. At the same time, use the temperature sensor integrated in the probe array to collect real-time temperature data at the bottom sediment measurement depth. Based on the empirical formula for temperature-resistivity of water-saturated sediment, temperature compensation is performed on the resistance value after deducting the influence of the polarization effect, and it is corrected to the resistivity value at the standard temperature of 25°C. The temperature compensation formula is: where ρ 25 is the resistivity at 25°C (unit: ohm·m), is the resistivity (in ohm·m) at temperature T (in °C), and 0.02 is the empirical temperature coefficient (in °C-1). The actual resistivity value of each measurement point, after polarization effects and temperature compensation, is associated with the corresponding timestamp, probe array ID, and measurement depth to form a set of calibrated resistivity values. Receive the set of calibrated resistivity values. For each probe array's resistivity data set at a specific time point (the set contains the calibrated resistivity values ​​of 12 electrodes at different depths on the probe array), the resistivity values ​​of these discrete depth points are interpolated using the cubic spline interpolation algorithm. Using the sediment surface as the depth zero point, interpolation is performed vertically downward 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 calculation results output the resistivity value at a fixed depth interval of 1 cm. The interpolation results of each probe array at each time point are saved to form a depth distribution curve dataset containing a timestamp, probe array ID, and a resistivity-depth sequence at the probe location at 1 cm intervals. The depth distribution curve data of all probe arrays at different time points are received. These data are organized according to the spatial location (latitude and longitude coordinates) and time of the probe arrays. To construct the three-dimensional resistivity distribution of the sediment area, the discrete probe array data are spatially interpolated using the Kriging interpolation method. The entire sediment monitoring area is divided into a grid with a horizontal resolution of 5 m × 5 m. For each grid cell, the depth distribution curve data of the nearby probe arrays are combined with the Kriging interpolation algorithm to calculate the resistivity prediction value of each depth layer (at 1 cm intervals) within the grid cell. The three-dimensional resistivity distribution field data calculated at different time points (for example, every hour) are stored in a time series. Finally, a sediment resistivity layer spectrum is formed, which records the resistivity distribution state of each depth layer of the sediment (1 cm resolution) at different spatial positions (5 m × 5 m grid resolution) and its changing characteristics over time (1 hour resolution), forming a resistivity distribution data set containing information in three dimensions: time, spatial position and depth.

[0039] Preferably, step S2 includes:

[0040] Step S21: pre-processing the reservoir environmental data to obtain a multi-source environmental parameter time series table;

[0041] Step S22: performing a lag correlation analysis on the sediment resistivity layer spectrum and the multi-source environmental parameter time series table to obtain a factor-depth response matrix;

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

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

[0044] In this embodiment of the present invention, water quality parameters (e.g., pH, dissolved oxygen concentration, turbidity, and water temperature) are collected through the reservoir's existing automated monitoring stations. Meteorological data (e.g., rainfall, air temperature, and air pressure) are collected using an automatic weather station. Hydrological data (e.g., water level using a water level gauge and inflow and outflow using a flow meter) are also collected. These sensors and equipment 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 greater than one hour, the data is aggregated into hourly values ​​using a mean aggregation method. If the raw data collection frequency is less than one hour, the data is interpolated into hourly values ​​using a linear interpolation method. This ensures that the time series of all environmental parameters are strictly aligned with the time points of the sediment resistivity layer spectra, with the hour as the minimum time unit. Missing data are then filled. For missing data at less than three consecutive time points, the average of the values ​​at the adjacent time points is used to fill the gap. Missing data at three or more consecutive time points are marked as invalid and excluded from subsequent analysis. Outliers were then removed by applying a median filter with a window length of five time points to the time series data for each environmental parameter. Values ​​that deviated from the median in the window by more than three standard deviations were marked as outliers and removed. Finally, the processed environmental parameter data were normalized using the Z-score method to transform the time series of each parameter into a distribution with a mean of 0 and a standard deviation of 1. The normalization formula is: Z = (X - μ) / σ, where X is the original value, μ is the mean of the parameter, and σ is the standard deviation of the parameter. The normalized time series data for each environmental factor were integrated into a table, recording all environmental parameter values ​​for each hour, to form a multi-source environmental parameter time series table. Sediment resistivity layer spectra (containing resistivity data at 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 layer spectra, the resistivity time series for each probe array position (or after spatial averaging) at each depth layer (e.g., every 1 cm depth) were first extracted. Then, the lagged correlation coefficient 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 the sediment resistivity at each specific depth layer is calculated. The lag time window is set from 0 to 72 hours, that is, the Pearson correlation coefficient between the environmental parameter time series and the resistivity time series is calculated at lags of 0 hour, 1 hour, ..., 72 hours. The correlation coefficient calculation formula is: Where r is the correlation coefficient, xi is the value of the environmental parameter at time point i, is the mean of the environmental parameter time series, y i+ τ is the value of resistivity at time point i + τ (lag time), is the mean of the resistivity time series, and ∑ represents the sum. For environmental parameters with significant seasonal variations (such as water temperature and water level) or resistivity changes, the STL decomposition method was used to decompose the time series into seasonal terms, trend terms, and residual terms. Lagged correlation analysis was then performed on the residual terms to eliminate seasonal interference. A sliding window method (with a window length set to 7 days) was used to calculate the stability of the lagged correlation coefficient over different time periods. That is, the standard deviation of the correlation coefficient within each window was calculated. The smaller the standard deviation, the more stable the correlation. Factor-depth pairs with an absolute value of the correlation coefficient greater than 0.5 and a sliding window standard deviation less than 0.1 throughout the monitoring period were 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 inverse of the standard deviation) were recorded to form a factor-depth response matrix, with rows representing environmental factors and columns representing sediment depth layers. Based on the factor-depth response matrix, a correlation network diagram of environmental factors was constructed. The nodes of the network diagram represent each environmental factor, and the lines between the nodes represent the correlation strength between the factors (the average correlation between the factors in the factor-depth response matrix or the correlation at a specific depth can be used). Based on this network diagram, the partial correlation analysis method is used to calculate the partial correlation coefficient between any two environmental factors, while controlling the influence of all other environmental factors in the network. Partial correlation analysis can reveal the direct correlation strength between factors and exclude indirect effects. The calculation results form a factor direct influence matrix. Analyze the factor direct influence matrix, identify factor combinations with significant positive partial correlation (synergistic enhancement effect) or negative partial correlation (antagonistic weakening effect), and classify these combinations and their interaction types to obtain a factor interaction pattern table. For example, high rainfall and high inflow have a synergistic effect. According to the factor interaction pattern table and the lag correlation analysis results in step S22 (especially the lag time information), the time lag characteristics of the combined effect of different factor combinations on the sediment resistivity are analyzed to obtain a time lag interaction matrix, which records the key factor combinations and their optimal time lags for affecting the sediment resistivity. Combining historical monitoring data with known sediment instability events (e.g., from historical disaster records), the numerical characteristics of the key factor combinations identified in the factor interaction pattern table before these events occurred are analyzed. Through statistical analysis or machine learning methods (e.g., support vector machine classifier), the critical threshold intervals or proportional relationships of these factor combinations when causing significant resistivity changes or signs of instability in the sediment are determined to obtain 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 factor-depth response matrix, the performance differences of these interaction patterns and conditional thresholds at different depth layers are analyzed to identify the characteristics of specific depth layers that are more sensitive to certain factor interaction combinations, and to obtain 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 vertical distribution of these interactions in the sediment, thereby revealing the complex mechanisms by which environmental changes influence sediment stability. Using the factor interaction effect map (depicting the interaction patterns, critical conditions, and depth distribution between factors) and the factor-depth response matrix (describing the impact and hysteresis of individual factors at different depths) as input, combined with the time points of sediment instability events marked in historical monitoring data, the changing characteristics of environmental factors and the changes in the sediment resistivity layer spectra in the period before these events were retrospectively analyzed. A discriminant model was trained using a decision tree learning method or logistic regression model, 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 signs of sediment instability (e.g., a rapid decrease in resistivity, a decrease in safety factor in stress field calculations, etc., defined based on historical experience or preset thresholds) as output labels. Discriminative rules are extracted from the trained decision tree model. These rules describe the logic for determining sediment stability when environmental factors reach specific numerical combinations or exceed specific thresholds. Furthermore, the model extracts the importance weights of each environmental factor (both single and combined) in determining sediment stability. The weights reflect the degree of influence of each factor on sediment stability. Based on the discriminative rules and historical data, the critical threshold intervals for 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 drawdown rate and pore water pressure change rate) are extracted. Combined with the factor-depth response matrix, the strength and manifestation of these critical conditions at different depths are determined. Ultimately, a sediment stability influencing factor spectrum is generated. This spectrum, in a structured format, records the influence weight of each key environmental factor (or factor combination) on sediment stability, the critical threshold interval for triggering instability (including numerical range or change rate), and the action characteristics or sensitivity differences of each factor (or combination) at different sediment depths (e.g., shallow, middle, and deep). This provides a quantitative influencing mechanism model for subsequent reconstruction of the sediment internal stress field and prediction of catastrophic criticality.

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

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

[0047] Perform partial correlation analysis on the factor correlation network diagram to obtain the factor direct influence matrix;

[0048] Classify and identify the interaction patterns of the factor direct influence matrix to obtain the factor interaction pattern table;

[0049] The time-lag interaction effect is analyzed based on the factor interaction pattern table to obtain the time-lag interaction matrix;

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

[0051] Perform depth difference analysis on the factor-depth response matrix according to the condition threshold table to obtain the depth interaction distribution map;

[0052] Draw a factor interaction effect diagram based on the depth interaction distribution diagram, factor interaction pattern table, and conditional threshold table.

[0053] In an embodiment of the present invention, a multi-source environmental parameter time series table is obtained, which contains standardized values ​​of environmental factors such as water quality parameters (pH value, dissolved oxygen, turbidity, water temperature), meteorological data (rainfall, temperature, air pressure) and hydrological data (water level, inflow, outflow) in 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 calculation of the correlation coefficient covers the entire monitoring period. Construct an undirected weighted graph, in which the nodes in the graph represent each environmental factor, and the lines between the nodes represent the correlation between the two factors. The weight of the line is set to the absolute value of the Pearson correlation coefficient of the corresponding two factor time series. This graph is drawn using a visualization tool, in which the thickness or color depth of the line reflects the size of the absolute value of the correlation coefficient, thereby intuitively showing the covariation relationship between environmental factors. The input is a multi-source environmental parameter time series table. For any pair of environmental factors F i and F j , calculate the value of the product under the control of all other environmental factors (i.e., excluding other factors affecting F i and F j After the linear influence of the common influence part, F i and F j The partial correlation coefficients between the two factors are obtained by calculating the ratio of the corresponding elements in the precision matrix. The partial correlation coefficients between all environmental factors are calculated and 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 is ​​the factor F. i and F jThe partial correlation coefficients between the two factors are shown in Table 1. This matrix reflects the direct correlation strength between environmental factors after removing indirect effects. The partial correlation coefficient values ​​in the matrix are analyzed. A threshold for the partial correlation coefficient is set. For example, if the absolute value is greater than 0.7, it is considered that there is a significant direct interaction relationship. If the partial correlation coefficient is greater than the threshold (and is positive), it is judged that there is a synergistic enhancement (Synergistic) effect between the factor pair, that is, their changes are 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 is negative), it is judged that there is an antagonistic weakening (Antagonistic) effect, that is, their changes are in opposite directions, and each will weaken the other's impact on sediment stability. Based on these judgment results, factor combinations with significant synergistic or antagonistic effects are identified (which can be pairwise combinations or multi-factor combinations), 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 to form a factor interaction pattern table. For example, a significant positive partial correlation between rainfall and inflow was identified, labeled a synergistic enhancement pattern; a significant positive partial correlation between the rate of change of water level and the rate of change of pore water pressure was also identified, labeled a synergistic enhancement pattern. For each interaction pattern identified in the factor interaction pattern table (e.g., the combination of factor A and factor B in a synergistic enhancement pattern), a composite time series representing the effect of the combination was first constructed based on its pattern type. For example, for synergistic enhancement, the standardized time series of factors A and B can be added together to form a composite time series S = standardized(A) + standardized(B). 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 layer spectrum was then 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 in each interaction pattern, the depth layer it affects, and the corresponding optimal time lag were recorded in a table to form a time lag interaction matrix. For example, the time-lag interaction matrix shows that the synergistic effect of rainfall and inflow is strongest at a depth of 0-30 cm in the sediment layer, with an optimal time lag of 6 hours; the synergistic effect of the rate of change of water level and the rate of change of pore water pressure is strongest at a depth of 50-80 cm in the sediment layer, with an optimal time lag of 12 hours. For each factor combination, impact depth layer, and optimal time lag τ recorded in the time-lag interaction matrix, the numerical characteristics of these factor combinations at time τ before the historical instability event are retrospectively analyzed. For example, if a historical landslide occurred at time T and the time-lag interaction matrix shows that the factor combination (A, B) has a significant effect at depth D and the optimal time lag is τ, the values ​​of factors A and B at time T-τ are analyzed.Through statistical analysis or rule mining, characteristic numerical ranges or relative relationships of these factor combinations are identified at time τ before an instability event occurs. For a synergistic enhancement model, the instability risk is identified to increase significantly when both factor A and factor B exceed a certain threshold at time T-τ (e.g., A > A0 and B > B0). For an antagonistic weakening model, the risk increases when one factor changes significantly while the other fails to provide sufficient counterbalance (e.g., C increases significantly while D does not decrease significantly). The critical numerical conditions or thresholds for each identified factor combination at a specific depth layer and time lag are recorded in a table to form a conditional threshold table. For example, the conditional threshold table may specify that the shallow layer risk increases when rainfall (6 hours ago) is greater than 20 mm / hour and the inflow change rate (6 hours ago) is greater than 10 cubic meters / second / minute; and the middle layer risk increases when the water level change rate (12 hours ago) is greater than 0.1 meters / hour and the pore water pressure (12 hours ago) is greater than 0.05 MPa. Analyze the distribution of critical conditions recorded in the condition threshold table across different depths. For example, some critical conditions are only effective within specific depth ranges. Furthermore, by combining the influence strengths and optimal time lags of individual factors at different depths in the factor-depth response matrix or the composite factors in the time-lag interaction matrix, identify differences in the vertical significance or sensitivity of different interaction patterns and their critical conditions within the sediment. For example, the critical condition for a particular interaction pattern may be lower at shallow depths but require higher values ​​to trigger at depth; or the influence strength of a particular interaction pattern may peak at a specific depth range. These depth-dependent 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 patterns, corresponding critical conditions, and impact strength assessments for different depth ranges (e.g., 0-30 cm, 30-80 cm, 80-180 cm); or a graphical representation showing the activity or importance of different interaction patterns along the depth axis. Combining all information from the depth interaction distribution map, the factor interaction pattern table, and the condition threshold table allows for the construction of a complete and interpretable model or map, known as a factor interaction effect map. This map clearly shows: which environmental factors interact with each other (from the factor interaction pattern table); the type of interaction (synergistic or antagonistic); the critical conditions or threshold combinations that trigger signs of instability (from the condition threshold table); and the manifestation characteristics and impact range of these interactions and critical conditions at different depths in the sediment (from the depth interaction distribution map). For example, the factor interaction effect map can be a flow chart or network diagram, with nodes representing environmental factors and edges representing interactions. The interaction type and critical conditions are labeled on the edges, and different edges or nodes are assigned different colors or styles to indicate their significance at different depths. This map provides a summary and quantitative description of how environmental factors complexly affect sediment stability.

[0054] Preferably, constructing the sediment layered structure model in step S3 includes:

[0055] The sediment physical parameters are inverted based on the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum to obtain the sediment physical property profile;

[0056] According to the resistivity layer spectrum of the sediment, the resistivity layer peak value of the sediment physical property profile is detected to obtain the layer interface candidate point set;

[0057] Perform characteristic clustering and stratification on the candidate point set of the layer interface to obtain a physical characteristic stratification scheme;

[0058] Identify the distribution map of the interlayer transition zones of the sediment physical property profile according to the physical property stratification scheme;

[0059] Identify the key weak layer map of the distribution map of the transition zone between layers;

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

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

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

[0063] In the embodiment of the present invention, the input is a sediment resistivity layer spectrum (containing resistivity data at each spatial position, depth, and time) and a sediment stability influencing factor spectrum (containing the weights and critical conditions of the influence of environmental factors on sediment stability). For each spatial position and time point in the sediment resistivity layer spectrum, its vertical resistivity profile is extracted. Based on the improved Archie formula, the porosity parameter is inverted using the resistivity data. The improved Archie formula adopts ρ = a × ρ_w × φ -m , where ρ is the resistivity of the sediment (unit: ohm·m), ρ_w is the resistivity of the pore water (estimated from the conductivity or total dissolved solids (TDS) in the water quality parameters, unit: ohm·m), φ is the porosity of the sediment (dimensionless), a is the lithologic coefficient (range: 0.5 to 2.5, with a value of 1 for pure sediment), and m is the cementation index (range: 1.3 to 2.5, with a value of 2 for loose sediment). The values ​​of a and m are set according to the sediment type (preliminarily determined by historical sedimentary data or the morphology of the resistivity profile). 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 solid particle density (range: 2.65 to 2.75 g / cm3 for sediment). Dry density ρ_d = (1-φ) × ρ_s. Total density ρ_bulk=ρ_d+φ×ρ_w. For cohesion c and internal friction angle φ i, combined with the sediment type and moisture content to estimate. For example, for silty clay, cohesion c is negatively correlated with moisture content w, and the internal friction angle φ is negatively correlated with moisture content w. i The estimated cohesion and internal friction angle are negatively correlated with the water content w. This relationship is represented by a piecewise linear or exponential function, whose parameters are calibrated based on historical field test data (e.g., cross-plate shear tests). Furthermore, the estimated cohesion and internal friction angle are corrected based on the influence weights and critical conditions of environmental factors (e.g., pore water pressure and water level change rate) on sediment strength as described in the sediment stability influencing factor spectrum. The estimated strength parameters are particularly reduced when the environmental factors approach critical values. The values ​​of density, water content, porosity, cohesion, and internal friction angle for 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 layer spectrum. For each spatial location and time point, the first-order derivative of the resistivity-depth profile curve with respect to depth, dρ / dz, is calculated. The derivative is calculated using the central difference method: (ρ(z+Δz)-ρ(z-Δz)) / (2Δz), where Δz is 1 cm. Before calculation, the original resistivity profile is smoothed using a Savitzky-Golay filter with a window length of 5 cm to reduce the influence of noise. The local maximum and minimum points in the dρ / dz curve (i.e., the zero-crossing point of the second-order derivative) are detected. The depths corresponding to these local extreme points are regarded as the locations where the resistivity change rate is the largest, usually corresponding to the layer interface where the properties of the bottom mud material or the water content change significantly. A threshold for the absolute value of the derivative is set, for example, |dρ / dz|>0.05 ohm·m / cm, and only the extreme points whose absolute value of the derivative exceeds this threshold are retained to exclude small fluctuations. The depths corresponding to these local extreme points exceeding the threshold are recorded as the candidate point set for the layer interface at that spatial position and time point.

[0064] The input is a set of candidate points of the layer interface and the profile data of the physical properties of the sediment. For each spatial location and time point, the vertical depth range of the location is divided into several preliminary layer segments using the candidate points of the layer interface. For each layer 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 of the layer segment. The feature vectors of all spatial locations, all time points, and all layers are collected. These feature vectors are clustered using the K-means clustering algorithm. The number of clusters K is set to 3, representing the main sediment types such as silt layer, silty silt layer, and silt layer. This number is set with reference to historical geological exploration data. The clustering results assign each layer segment to a specific physical property category (i.e., cluster cluster). Layer segments of the same category and adjacent in the vertical direction are merged to form a layer with relatively uniform physical properties. For each spatial location and time point, the vertical layer structure, including the thickness, burial depth, and physical property category of each layer, is determined based on the clustering results to form a physical property stratification scheme. The input is 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) identified in the physical property stratification scheme, a 10 cm depth range is extended above and below the interface as an inspection window. Within this window, the gradient of the sediment physical property profile data (e.g., water content or density) is analyzed. If the absolute value of the gradient of any physical parameter (e.g., water content) within the window exceeds 20% of the difference between the average parameters of the two major 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 / cm3 / cm3, the area is considered an interlayer transition zone. The starting and ending depths of the transition zone are determined, i.e., the depth range that meets the gradient threshold condition. The spatial position, time point, depth range of the identified interlayer transition zone and the corresponding physical property change characteristics are recorded to form an interlayer transition zone distribution map. The input is 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 φ in the interlayer transition zone are extracted from the sediment physical property profile data. i According to the preset strength threshold (for example, cohesion c<3kPa or internal friction angle φ i<15 degrees), marking transition zones with lower intensity. Then, referring to the spectrum of factors influencing sediment stability, depth intervals particularly sensitive to changes in pore water pressure or stress were identified. The lower intensity transition zones were overlaid with depth intervals sensitive to environmental changes. If an interlayer transition zone has both low intensity and high sensitivity to environmental changes (e.g., an absolute correlation with pore water pressure >0.7), it was marked as a critical weak layer. The spatial location, time point, depth range, and weakness characteristics (low intensity, high sensitivity) of each identified critical weak layer were recorded to form a critical weak layer map. The inputs were a physical property stratification scheme and a digital elevation model (DEM) of the reservoir bottom. The reservoir bottom DEM was acquired using a high-precision single-beam or multi-beam bathymetry 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 vertical depth of the layer interface (relative to the sediment surface) at each monitoring location (discrete point). It was assumed that layer interfaces within the same physical property category are continuous and smoothly varying horizontally and are roughly parallel to the sediment bottom topography. For each grid point (x, y) on the DEM, first obtain the corresponding sediment bottom elevation Z_bottom(x, y). Then, find the three monitoring locations closest to the grid point and obtain the physical property stratification scheme of these three locations. Based on the relative depth information of the layer interface of these three locations and their horizontal distances to the (x, y) point, use the inverse distance weighted method or kriging interpolation method to estimate the depth d of the layer interface of each physical property category at the (x, y) point relative to the sediment surface. i (x,y). The position of the layer interface at the 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 boundaries (i.e., the layer interface elevation distribution) of each physical property category layer within the entire monitoring area. These boundaries together constitute a three-dimensional layered terrain model of the sediment, which describes the layered geometry of the sediment in space. The input is a three-dimensional layered terrain model and a sediment physical property profile dataset. The three-dimensional layered terrain 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, and internal friction angle) of all layers belonging to that category are collected from the sediment physical property profile dataset. 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 terrain model. The interpolation calculation estimates the corresponding physical parameter value for each cell (e.g., finite element grid cell) in the three-dimensional layered terrain model. The interpolation process takes into account the spatial correlation of the parameters. The resulting layered physical parameter field is a three-dimensional dataset in which each spatial location (or grid cell) is assigned a set of spatially interpolated physical parameter values, such as density, water content, porosity, cohesion, and internal friction angle, consistent with its physical property category. The inputs are the layered physical parameter field and a critical weak zone map. The layered physical parameter field already provides the three-dimensional physical property distribution of the bulk of the sediment. Information from the critical weak zone map is superimposed on this three-dimensional model. The critical weak zone map identifies weak areas within specific spatial locations and depth ranges. When generating the layered sediment model for subsequent finite element analysis, the volume elements corresponding to these critical weak zones are specifically marked. Based on their degree of weakness (e.g., based on the assessment of weakness characteristics from the critical weak zone map), these areas can be assigned modified, lower strength parameter values ​​(cohesion, internal friction angle), or locally refined during finite element meshing. The resulting sediment layered structure model is a digital representation that includes the three-dimensional geometry of the sediment, its internal layered structure, and the spatially varying physical parameters of each layer and location. It also clearly identifies and highlights key weak layers, providing precise geometric and material property inputs for subsequent stress field calculations.

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

[0066] The stress-strain relationship of the sediment layered structure model is calculated to obtain the sediment mechanical response parameter set;

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

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

[0069] In the embodiment of the present invention, the input is the sediment layer structure model, which provides the geometric layering information of the sediment in space and the physical parameters of each layer at each position (density ρ_bulk, porosity φ, water content w, cohesion c, internal friction angle φ i Based on these physical parameters, the stress-strain constitutive relationship of the sediment material is defined. An elastoplastic constitutive model, such as one based on the Mohr-Coulomb yield criterion, is used. This model follows Hooke's law in the elastic phase 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 from the sediment density ρ_bulk and porosity φ using an empirical formula, for example, E = k × (ρ_bulk) 2 / φ p , where k and p are empirical coefficients whose values ​​are determined based on the sediment type (determined by the physical property stratification scheme) and historical field test data. Poisson's ratio ν is usually taken as a value close to 0.5 for saturated soft soils, such as 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 layer structure model. i ) together constitute the description of the mechanical response characteristics of the sediment at different locations and under different conditions. By associating these parameter sets with the spatial positions in the sediment hierarchical structure model, the sediment mechanical response parameter set is formed, which is a spatially distributed parameter field.

[0070] The input is the sediment hierarchical structure model and the sediment mechanical response parameter set. Using a three-dimensional finite element pre-processing tool, discretization is performed using a structured or unstructured mesh generation algorithm based on the geometric boundaries of the sediment hierarchical structure model. Three-dimensional solid elements (e.g., hexahedral elements or tetrahedral elements) are generated. Horizontally, the side length of the mesh element is set between 1 and 5 meters, with the specific value determined based on the reservoir area size and the required computational accuracy. Vertically, the height of the mesh element is set between 5 and 20 centimeters to accurately capture stress variations in the vertical direction. Based on the key weak layers identified in the sediment hierarchical structure model (from the substep of constructing the sediment hierarchical structure model in S3), the mesh in these areas is locally refined to reduce the mesh size. For example, the vertical mesh size can be reduced to 2 centimeters. After mesh generation is complete, all generated finite element mesh elements are traversed. For each mesh element, its spatial position (e.g., the coordinates of the cell center) is determined. Based on the layer and physical property category to which this position belongs in the sediment hierarchical structure model, the corresponding physical parameters (E, ν, c, φ) in the sediment mechanical response parameter set are searched. i These parameters are assigned to the mesh element as material properties. The resulting sediment finite element mesh model is a digital representation of the node coordinates, element connectivity, and the material properties (mechanical response parameters) of each element.

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

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

[0073] The dynamic permeability coefficient field is obtained by dynamically correcting the layered permeability coefficient of the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum.

[0074] Perform multi-source water pressure superposition calculation based on the dynamic permeability coefficient field to obtain a pressure component analysis table;

[0075] Identify the interface pressure jump according to the pressure component analysis table and obtain the pressure jump position map;

[0076] Analyze the pressure wave velocity distribution diagram based on the pressure component analysis table;

[0077] The critical instability pressure is predicted based on the pressure wave velocity distribution map and the pressure jump position map to obtain the pore pressure risk map;

[0078] Generate a pore water pressure distribution map based on the initial stress distribution map and the pore pressure risk map;

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

[0080] In the embodiment of the present invention, the input is a sediment resistivity layer spectrum (including resistivity data at each spatial position, depth and time) and a sediment stability influencing factor spectrum (including the influence weights and critical conditions of environmental factors on sediment stability). First, based on the resistivity values ​​in the sediment resistivity layer spectrum, the porosity φ of the sediment is inverted and calculated using the empirical relationship based on the Archie formula, for example, φ = (a×ρ_w / ρ)^(1 / m), where ρ is the resistivity, ρ_w is the pore water resistivity, and a and m are empirical coefficients, whose values ​​are set according to the sediment type (determined by the sub-step of constructing the sediment hierarchical structure model in step S3). Then, based on the porosity φ and sediment type obtained by inversion, the Kozeny-Carman equation is used to estimate the initial permeability coefficient k0 of the sediment: k0 = C×φ 3 / (1-φ) 2, where C is a constant related to the particle shape and size distribution, and its value is determined by the type of sediment. This k0 is a permeability coefficient based on static physical properties. Next, referring to the spectrum of factors affecting sediment stability, identify environmental factors that have a significant impact on pore water pressure changes (e.g., water level change rate, rainfall intensity), their critical thresholds, and influence weights. When these environmental factors approach or exceed the critical threshold, the permeability of the sediment will change dynamically (e.g., rapid water level drop will cause surface cracks to increase permeability, and sustained high water pressure will cause particle migration to block pores and reduce permeability). Based on the deviation of the current value of the environmental factor relative to its critical threshold and its influence weight, a dynamic correction factor F_dynamic is calculated. 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 more it exceeds the threshold, the larger the function value). The dynamic permeability coefficient k = k0 × F_dynamic. The dynamic permeability coefficient k calculated at each spatial position, 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 is the dynamic permeability coefficient field and the multi-source environmental parameter time series table (including water level changes, rainfall and other data). Using the three-dimensional transient seepage finite element model, the sediment volume is discretized into a grid compatible with the stress field calculation. The dynamic permeability coefficient field is input into the seepage model as a material property. According to the multi-source environmental parameter time series table, time-varying boundary conditions are imposed on the boundaries of the seepage model: a time-varying head boundary related to the reservoir water level height is imposed on the sediment surface (reservoir bottom) (head h = water level elevation); a flow boundary related to rainfall intensity is imposed on the infiltration boundary that may exist (such as the rainfall area); and a zero flow boundary is imposed on the impermeable boundary (such as the bedrock interface). Solve the transient seepage equation. Where ρ_w is the density of water, S_s is the water storage rate (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 and sink term (such as rainfall infiltration). Solve to obtain the time-varying hydraulic head h(x,y,z,t) at each point inside the sediment. Convert the hydraulic head to pore water pressure u(x,y,z,t) = ρ_w × g × h(x,y,z,t), where g is the acceleration of gravity. Analyze the calculated pore water pressure field. It can be decomposed into a hydrostatic pressure component (pressure caused only by the current water level) and an excess hydrostatic pressure component (additional pressure caused by transient effects such as seepage, water level changes, and rainfall). Record the total pore water pressure value and (optionally) its decomposed components at each spatial location, depth, and time point in a table to form 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 position and time point, the pore water pressure distribution curve u(z) is analyzed along the vertical depth direction. The gradient of pore water pressure with respect to depth is calculated. Identification There is a significant mutation in the curve (i.e. the second derivative The depth position where there is a local extreme value or a change in sign. These locations usually correspond to the interface of the sediment layer or the place where the permeability coefficient changes significantly. Set the threshold for the sudden change of pressure gradient, such as The rate of change of the pressure jump exceeds a certain percentage threshold, or the absolute value of the second-order derivative exceeds a certain threshold. The depth locations with significant pressure jumps are recorded along with their corresponding spatial coordinates and timestamps 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). Typical transient events occurring in the multi-source environmental parameter time series are identified, such as a rapid drop in water level or a significant rainfall event. In the pressure component analysis table, the propagation of pore-water pressure disturbances caused by these events within the sediment is tracked. For a specific spatial location, the changes in the pore-water pressure time series at different depths are analyzed to identify the times when the pressure peaks or wavefronts reach different depths. The average velocity of the pressure wave propagating from one depth to another is calculated as Δz / Δt, where Δz is the depth difference and Δt is the time required for the pressure wave to propagate. This process is repeated to calculate the pressure wave velocity at different spatial locations and depth intervals (e.g., shallow, mid-layer, deep) under the influence of different environmental events. The calculated pressure wave velocity and its corresponding spatial location, depth interval, and time are recorded to form a pressure wave velocity distribution map. This map reflects the propagation speed of pore pressure changes within the sediment and is closely related to the permeability and compressibility of the sediment.

[0084] The inputs are a pressure velocity distribution map, a pressure jump location map, and a spectrum of factors influencing sediment stability. The pressure velocity distribution map is analyzed to identify areas of abnormal pressure velocity (for example, excessively fast velocities indicate a dominant seepage channel or a high permeability region, while excessively slow velocities indicate reduced permeability or stress concentration). The pressure jump location map is analyzed to identify interfaces with significant and persistent pressure jumps. This information, combined with the information in the sediment stability factor spectrum, indicates that instability can occur when pore water pressure reaches a critical value or when the pore water pressure ratio reaches a critical value (for example, u / σ_z approaches or exceeds 1). The current pore water pressure u(x,y,z,t) is compared with the total stress σ_z(x,y,z) at that location (derived from the initial stress distribution map or taking into account subsequent stress changes) to calculate the pore water pressure ratio r_u = u / σ_z. Areas where r_u approaches or exceeds the critical threshold are identified. Areas of high r_u, interfaces with significant pressure jumps, and areas of abnormal pressure velocity are superimposed for analysis. If an area simultaneously meets multiple abnormal conditions (high r_u, pressure jump, and abnormal velocity), the pore water pressure in that area is considered to be close to a critical instability state and has a high risk level. Risk scores or grades are assigned to each sediment area based on the r_u value, the pressure jump amplitude, and the degree of velocity anomaly. The spatial distribution of risk scores or grades is represented to form a pore pressure risk map, which visually illustrates areas 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 previous step already calculated the pore-water pressure field u(x,y,z,t) (in the Pressure Component Analysis table), the description of 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 could mean using the calculated pore-water pressure field u(x,y,z) (at the current point in time) and verifying, revising, or highlighting the pore-water pressure values ​​in high-risk areas based on the information in the pore pressure risk map. For example, if the calculated pore-water pressure indicates a high-risk area in the risk map, that area will be specifically marked or colored in the final pore-water pressure distribution map. Alternatively, this step aims to ultimately output a clear, easy-to-understand spatial distribution map of sediment pore-water pressure at the current moment, u(x,y,z), supplemented with risk information. Considering that the subsequent steps require accurate pore water pressure values ​​to calculate effective stress, the most reasonable explanation is that this step outputs the accurate pore water pressure spatial distribution map u(x,y,z) at the current moment, which has been verified or confirmed by the previous analysis (including risk assessment).

[0086] The input is the pore water pressure distribution map (including the pore water pressure u(x,y,z) at each point in the sediment at the current moment) and the initial stress distribution map (including the total stress σ(x,y,z) at each point in the sediment). According to the Terzaghi effective stress principle, for each point in the sediment, the effective stress tensor σ'(x,y,z) is calculated. The normal stress component (σ x ,σ γ The effective stress of σ2) is obtained by subtracting the pore water pressure from the total normal stress: x =σ x -u,σ' γ =σ γ -u, σ'2=σ2-u. Shear stress component (τ xγ ,τ γ2 ,τ 2x ) is not directly affected by pore water pressure, and its effective stress is equal to the total shear stress: τ' xγ =τ xγ , τ' γ2 =τ γ2 , τ' 2x =τ 2x In the calculation, all stress components and pore water pressures must use a consistent sign convention (for example, compressive stress is positive or negative). The effective stress distribution diagram is formed by storing the information of the stress and its spatial position. This diagram reflects the actual stress state of 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] Compute a multi-scale gradient tensor set for effective stress distribution maps;

[0089] Compute stress curl profiles based on a multi-scale gradient tensor set;

[0090] Calculate the critical gradient ratio map based on the multi-scale gradient tensor set and the sediment layer structure model;

[0091] Perform gradient time evolution processing on the multi-scale gradient tensor set to obtain the gradient evolution trend graph;

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

[0093] In the embodiment of the present 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 unit in the sediment finite element mesh model Data that changes over space and time. For each time point, the spatial gradient of the effective stress distribution map is calculated. Gradient operator In the three-dimensional Cartesian coordinate system, it is expressed as Apply the gradient operator to each component of the effective stress tensor, for example, calculate the gradient of the effective normal stress σ'x These gradient components together constitute the information of the stress tensor gradient. The calculation uses numerical differentiation methods, such as central difference or finite difference method based on finite element mesh. In order to achieve "multi-scale" calculation, before calculating the gradient, Gaussian smoothing filters with different standard deviations σ_s are applied to the effective stress component field. For example, a three-dimensional Gaussian filter with standard deviations of 0.5 meters, 1 meter, and 2 meters is used to convolve the effective stress component field, and then the gradient is calculated on each smoothed field. The standard deviation σ_s represents the scale of smoothing. Smaller σ_s retains more details, and larger σ_s focuses on macro trends. The stress component gradient vectors calculated at different scales (for example, for σ'x, obtained at scale i) are combined. ) are collected to form a multi-scale gradient tensor set, which contains the stress tensor gradient information at each spatial position, time point, and different scales.

[0094] The input is a set of multi-scale gradient tensors. Stress curl describes the local rotation or vortex tendency of the stress field and is closely related to the shear deformation and potential failure of the 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 These partial derivative information can be extracted from the multi-scale gradient tensor set (e.g., yes y component of the stress curl vector. At each selected scale, the three components of the curl vector are calculated based on the data in the multiscale gradient tensor set. The stress curl vector calculated at each spatial location, time point, and different scales is stored.

[0095] Usually we focus on the modulus of the curl vector

[0096] As the main content of the stress curl distribution map, it shows the rotation intensity of the local stress field. This map reflects the unevenness and vortex degree of the local shear stress field inside the sediment.

[0097] The input is a multi-scale gradient tensor set and a sediment layer structure model (the model provides the cohesion c and internal friction angle φ at each location i The critical gradient ratio is a measure of the severity of the stress gradient relative to the material strength. For each point within the sediment, calculate its local shear strength. in is the effective normal stress at the point (obtained from the effective stress distribution map). Select a stress gradient modulus of a representative scale from the multi-scale gradient tensor set, for example, select the stress tensor gradient norm at a specific scale (for example, 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 tensor components) or the maximum shear stress gradient modulus Calculate the critical gradient ratio or High R values ​​indicate large stress gradients in areas of 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 identifies areas within the sediment where the rate of stress change is abnormally high compared to the local strength. These areas are potential locations for stress concentration and failure initiation.

[0098] The input is a set of multi-scale gradient tensors (containing gradient information at different time points). For each spatial position and each selected scale inside the sediment, the stress gradient tensor (or its modulus, component) is analyzed over time. The rate of change of the stress gradient over time, i.e., the time derivative, is calculated. or The time derivative is calculated using the time series difference method, for example, the difference between the gradient values ​​of two time points is divided by the time interval. In order to smooth the time series noise, the gradient time series can be processed with a sliding average of a window length of 3 time points before calculating the time derivative. Identify the areas where the gradient increases rapidly over time. The time derivative of the gradient calculated at each spatial location and scale, or its indicated evolutionary trend (e.g., rapid increase, stability, rapid decrease), is stored to form a gradient evolution trend graph. This graph reveals the dynamic characteristics of the stress gradient within the sediment over time. A rapidly increasing gradient trend is a key precursor to impending instability.

[0099] The input is a critical gradient ratio map, a stress curl distribution map (focusing on the curl modulus), and a gradient evolution trend map. The information from these three maps is comprehensively analyzed to identify areas with abnormal stress gradients. The threshold conditions for judging abnormalities are set. For example, an area is marked as abnormal if it meets the following conditions at the same time: ① The critical gradient ratio R exceeds the threshold R_threshold (for example, R>0.8); ② The stress curl modulus Exceeds the threshold Curl_threshold (for example, Pascal / meter); ③ The gradient evolution trend indicates that the gradient is increasing rapidly, that is, Exceeding the threshold Rate_threshold (for example, Pascal / meter / hour). These thresholds are determined based on a historical disaster precursor feature library (concepts from step S4.1) or engineering experience. Identify continuous areas that meet abnormal conditions in three-dimensional space. According to the degree of satisfaction of the conditions, the size of the area, and the numerical value of the abnormal value, the abnormal area is graded (for example, first-level abnormality, second-level abnormality) or assigned an abnormal index. The spatial position, shape, abnormal level or index of the identified abnormal area is visualized to form a stress gradient abnormal area map. This map clearly points out the areas where the stress state inside the sediment is abnormal and is developing towards instability, and is a key basis for early warning.

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

[0101] Identify potential sliding surfaces based on the stress gradient anomaly area map and the intensity ratio distribution map to obtain a potential sliding surface distribution map;

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

[0103] The temporal and spatial evolution of the stress field is analyzed based on the effective stress distribution diagram, the sliding surface safety factor table and the stability influencing factor spectrum, and the stress field spectrum between the bottom mud layers is obtained.

[0104] In an embodiment of the present invention, the input is a stress gradient anomaly area 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, the two maps are superimposed in three-dimensional space. Grid cells or regions that simultaneously meet the following conditions are identified: their strength ratio is close to or exceeds a preset threshold (for example, the strength ratio is >0.7) and is located within the stress gradient anomaly area (for example, the anomaly index is >0.5). These areas indicate locations where the material strength is relatively low and the stress changes drastically, and are the starting points for potential damage initiation. Next, a connected domain analysis method based on graph theory is adopted. Grid cells that meet the above conditions are regarded as nodes of the graph. If two nodes are spatially adjacent and their connection direction is roughly consistent with the principal shear stress direction of the region (obtained from the effective stress distribution map) (for example, the angle is less than 30 degrees), a connecting edge is established between them. On the constructed connected graph, continuous connected areas of a certain scale (for example, containing more than 100 grid cells) are searched. These areas are regarded as potential sliding surfaces. To obtain a more accurate sliding surface geometry, surface fitting techniques (e.g., polynomial or spline surface fitting based on the least squares method) or minimum energy path search algorithms are applied to the nodes within the identified potential sliding surface area. The minimum energy path search searches for a path connecting a starting point (e.g., the point with the strongest stress gradient anomaly) and an end point (e.g., a free boundary) in an area with high stress gradient and high intensity ratio. The path with the lowest cumulative "energy" (e.g., a weight related to the intensity 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 grid cells constituting the surface or parameters of the fitted surface) and the characteristic values ​​associated with the surface (e.g., average intensity ratio, maximum stress gradient anomaly index) are recorded to form a potential sliding surface distribution map.

[0105] The input is a potential sliding surface distribution map (including the geometric information of the identified potential sliding surfaces). For each potential sliding surface identified in the potential sliding surface distribution map, the safety factor is calculated using the limit equilibrium method. The sliding body enclosed by the sliding surface is divided into a series of vertical or inclined blocks (for example, the blocks are divided into blocks at intervals of 1 meter along the horizontal direction). For each block, the average effective normal stress σ' at its bottom (i.e., the sliding surface) is obtained from the effective stress distribution map. and the average shear stress τ. From the sediment physical property profile (or layered physical parameter field), the effective cohesion c' and effective internal friction angle φ of the bottom material of the strip are obtained. i Calculate the anti-slip force of each bar where l i is the length of the bottom of the bar on the sliding surface. Calculate the sliding force T of each bar i , usually the deadweight of the bar W iThe component along the sliding surface may also require consideration of external loads (such as additional stress caused by changes in water pressure). The safety factor F_s of the entire sliding body is calculated according to a specific method of the limit equilibrium method (for example, the simplified Bishop method or the Spencer method). The safety factor formula of the simplified Bishop method is: F_s = [∑(c'×l i +(W i -u i ×l i )×tanφ i ) / m γi ] / ∑(W i ×sinα i ), where u i is the pore water pressure at the bottom of the strip (obtained from the pore water pressure distribution map), α i is the inclination angle of the sliding surface at the bottom of the block, m γi =cosα i +sinα i ×tanφ i / F_s (requires iterative solution). The Spencer method considers the balance of forces and moments between bars, resulting in more accurate calculations. Record the unique identifier of each potential sliding surface and the calculated safety factor in a table to form a sliding surface safety factor table.

[0106] The input is the effective stress distribution diagram (including the effective stress field σ'(x,y,z,t) that changes with time), the sliding surface safety factor table (including the safety factor F_s(t) of each potential sliding surface that changes with time), and the sediment stability influencing factor spectrum (including the influence mechanism of environmental factors on sediment stability). First, the time series data of the effective stress distribution diagram is analyzed. For key areas (for example, 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 sliding surface safety factor table. For each potential sliding surface, plot its F_s curve over time and calculate the time derivative of the safety factor dF_s / dt. Pay special attention to sliding 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 an unstable state. Combined with the spectrum of factors affecting sediment stability, the state of the current environmental factors (obtained from the multi-source environmental parameter time series table) is associated with the changes in the stress field and safety factor. For example, if the water level drops rapidly and dF_s / dt is negative, and the influence factor spectrum shows that the water level drop will reduce the pore water pressure and thus increase the effective stress (usually increase the safety factor), it is necessary to further analyze whether there are other antagonistic factors or abnormal responses of specific layers that cause the safety factor to decrease, or whether there are deviations in the current calculation model. Leveraging current trends in environmental factors and historical correlations between environmental factors and stress / safety factor variations, the effective stress trends at key locations and the safety factor trends at key potential sliding surfaces are predicted over the next period (e.g., the next 24 or 72 hours). This dynamic information—including the real-time effective stress distribution at each depth level of the sediment, the real-time pore water pressure distribution (obtained from the 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-dependent rates of change of key stress parameters and safety factors, and predicted future trends—is integrated and visualized to generate a stress field map for interlayer sediment layers. This map provides a comprehensive, dynamic representation of the internal mechanical state of the sediment and its characteristics as it evolves 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: Acquire historical disaster data of reservoir sediment; extract historical disaster features from the historical disaster data of reservoir sediment to obtain a disaster precursor feature library;

[0109] Step S42: performing stress field anomaly detection based on the stress field map between sediment layers to obtain a stress field anomaly area table;

[0110] Step S43: performing an instability probability assessment based on the stress field anomaly area table and the disaster precursor feature library 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 abnormal area table to obtain the sediment disaster critical index;

[0112] Step S45: generating and issuing warning information based on the sediment disaster critical index to obtain a reservoir sediment disaster warning report.

[0113] In an embodiment of the present invention, all records of sediment-related disaster events since the reservoir was built are collected, including but not limited to sediment landslides, mudflows, liquefaction, etc. These records can come from accident reports, monitoring data, on-site photos, videos, media reports and relevant research literature from the reservoir management department. For each disaster event, its time of occurrence, specific location, impact range, disaster type (e.g., shallow sliding, deep sliding, mudflow), triggering factors (e.g., heavy rain, sudden drop in water level, earthquake) and the losses caused are recorded in detail. If so, the monitoring data for a period of time before the disaster occurs are retrospectively queried, especially environmental data such as reservoir water level, rainfall, inflow, and (if any) data such as sediment surface displacement and pore water pressure recorded by the historical monitoring system. Using these data, combined with the method of stress field calculation and potential sliding surface analysis in step S3, the sediment stress field state and potential sliding surface characteristics before the historical disaster occur are simulated or inversely analyzed. Typical stress field precursor signals before a disaster occur are extracted, such as peak shear stress at key locations, rapid increases in stress gradients, anomalous increases in pore water pressure leading to a significant decrease in effective stress, a potential slip surface safety factor dropping to a critical value (e.g., between 1.1 and 1.3) and continuing to decline, and a rapid decrease in resistivity at specific depths. The extracted characteristic stress field precursor patterns for different types of sediment hazards (e.g., shallow slip precursors: high strength ratios concentrated in shallow layers and rapidly increasing shallow stress gradients; deep slip precursors: a continuous decrease in the safety factor of deep potential slip surfaces; mudflow precursors: a rapid increase in pore water pressure leading to near-zero effective stress) are structured and stored along with their corresponding environmental triggering conditions and occurrence time windows, forming a database of disaster precursor signatures. This database serves as a benchmark for subsequent anomaly detection and probability assessment.

[0114] The input is the current interlayer stress field map of the sediment (including real-time effective stress distribution, pore water pressure distribution, potential sliding surface location, safety factor, and other information). The current stress field map is analyzed based on the stress field precursor signal types and thresholds defined in the disaster precursor feature library. For example, the analysis checks whether there are areas within the sediment where the ratio of shear stress to shear strength exceeds 0.8; whether there are areas where the stress gradient (for example, 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 (as determined by the gradient evolution trend graph); whether there are areas where the pore water pressure ratio (u / σ_z) exceeds the threshold 0.9; and whether there are sliding surfaces in the potential sliding surface safety factor table with a safety factor less than 1.3 and a time derivative dF_s / dt less than -0.01 / hour. Any area or potential sliding surface that meets any of the abnormal conditions is marked as a stress field anomaly. Record the spatial location (e.g., center coordinates or sliding surface ID) of each anomalous region or potential sliding surface, the type of anomaly (e.g., intensity ratio anomaly, gradient anomaly, pore water pressure anomaly, safety factor anomaly), the degree of anomaly (e.g., magnitude exceeding a threshold), and the rate of change of the anomaly parameter. This information is compiled into a table to form a stress field anomaly region table. This table lists all anomalous signals currently present in the sediment that are similar to historical precursory patterns.

[0115] The input is a table of stress field anomaly regions and a database of disaster precursory features. For each anomaly region or potential slip surface in the stress field anomaly region table, its current stress field characteristics (anomaly type, degree, and rate of change) are compared with precursory patterns of different disaster types in the database. The similarity between the current anomaly characteristics 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 high-intensity ratio region with a rapidly increasing stress gradient is currently monitored in the shallow layer, it has a high similarity to the precursory pattern of a shallow slip disaster. Based on this similarity, the probability that the anomaly region or potential slip surface will eventually develop into an actual disaster is assessed. This assessment can be accomplished using Bayesian networks, support vector machines, or rule-based reasoning systems. For example, P(disaster|current anomaly characteristics) = P(current anomaly characteristics|precursorial pattern) × P(precursorial pattern) / P(current anomaly characteristics), where P(precursorial pattern) can be determined based on the frequency of historical disaster events. The size and spatial distribution of the anomaly region are also considered. If multiple interrelated abnormal areas appear at the same time, or the abnormal areas are located in key weak layers, the overall probability of instability increases. For potential sliding surfaces, the safety factor F_s is the core indicator. The probability of instability is nonlinearly negatively correlated with F_s and positively correlated with the absolute value of dF_s / dt. By analyzing historical data, a quantitative relationship curve between F_s and dF_s / dt and the probability of instability is established. The instability probability of each spatial position (or grid unit) is calculated to form a three-dimensional probability distribution field. The probability distribution field is visualized to form a regional instability risk map, in which different colors or grayscales are used to indicate the probability of instability in each area.

[0116] The input is a regional instability risk map and a stress field anomaly area table. The sediment catastrophe critical index is an indicator that comprehensively quantifies the current overall instability risk of the sediment. The calculation of this index comprehensively considers the following dimensions: ① Instability probability: extract the highest probability value or the average probability of the high-risk area from the regional instability risk map. ② Abnormality degree: extract the most serious abnormality type, the highest abnormality degree or the abnormality index from the stress field abnormality area table. ③ Spatial scope: calculate the proportion of the total volume or area of ​​high-risk areas (for example, instability probability > 0.6) to the entire monitoring area. ④ Time urgency: based on the rate of change of key abnormal parameters (such as safety factor, stress gradient), combined with the time required for these parameters from the appearance of the anomaly to the occurrence of the catastrophe in the historical precursor feature library, estimate the time window required to reach the critical state (for example, the safety factor drops to 1.05). The critical index calculation formula can be a weighted summation or multi-factor product. For example, the catastrophic critical index = W1 × P_max + W2 × A_max + W3 × S_ratio + W4 / T_window, where P_max is the maximum probability of instability, A_max is the maximum degree of anomaly, S_ratio is the area ratio of the high-risk area, and T_window is the expected time window to the critical state. W1, W2, W3, and W4 are weighting coefficients, whose values ​​are calibrated based on the analysis of historical disaster cases. Based on the numerical range of the critical index, graded warning standards are set. For example, an index < 20 is normal; 20 ≤ < 50 is caution; 50 ≤ < 80 is warning; 80 ≤ < 100 is danger; and 100 ≥ is emergency. The calculated value is the sediment catastrophic critical index, which reflects the overall risk level and urgency of sediment instability.

[0117] The input is the critical index of sediment disaster. According to the calculated critical index value, the current warning level (normal, attention, warning, danger, emergency) is determined. According to the current stress field abnormal area table and regional instability risk map, the location of the area with the highest risk, the potential disaster type (for example, if the abnormality is mainly concentrated in the shallow layer and has a high correlation with rainfall, it will be shallow sliding; if the deep safety factor continues to decline, it will be deep sliding) and the estimated impact range (for example, through the spatial range estimation of the high-risk area). According to the time window estimated in step S44, the time period in which the disaster is expected to occur is given. This information is integrated into a structured reservoir sediment disaster warning report. The report content at least includes: the current warning level, the warning release time, the spatial location description of the abnormal area (for example, the area name, the longitude and latitude range), the main abnormal characteristics (for example, the safety factor and change rate of the critical sliding surface, the highest intensity ratio area), the potential disaster type, the estimated impact range, the expected occurrence time window, and the recommended emergency response measures for the warning level and disaster type (for example, strengthening inspections, limiting the reservoir water level drawdown rate, activating emergency plans, notifying downstream, etc.). For warning, danger, and emergency level alerts, the system automatically sends alert notifications to pre-defined reservoir management department heads and on-duty personnel via various channels, including text messages, emails, and phone calls. Simultaneously, intuitive warning visualizations are generated, including marking risk areas on a reservoir map, displaying a time series graph of criticality indices to illustrate risk trends, and showing safety factor curves for key potential slip surfaces. This visualization information is published alongside the warning report, allowing managers to quickly understand the risk situation and take action. The final output is a reservoir sediment disaster warning report.

[0118] It is particularly important that the probability of instability be assessed as follows:

[0119] Identify the regional instability mode table based on the stress field abnormal area table;

[0120] According to the regional instability pattern table, the historical similarity calculation is performed on the stress field abnormal area table and the disaster precursor feature database to obtain a historical similar case matching table;

[0121] Based on the stress field abnormal area table and regional instability mode table, the shear sliding type instability safety assessment is carried out to obtain the sliding surface safety factor assessment table;

[0122] According to the stress field abnormal area table and regional instability mode table, the liquefaction type instability risk assessment is carried out to obtain the liquefaction risk assessment table;

[0123] The single-point instability probability is calculated based on the historical similar case matching table, sliding surface safety factor evaluation table, and liquefaction risk evaluation table to obtain the point instability probability table;

[0124] According to the point instability probability table, the spatial correlation analysis of the stress field map between the bottom mud layers is carried out to obtain the risk space correlation map;

[0125] The cascading effect is evaluated on the point instability probability table and the risk space association map to obtain the disaster diffusion risk map;

[0126] Conduct time window prediction on the disaster diffusion risk map to obtain the regional instability risk map;

[0127] In the embodiment of the present invention, according to the abnormal feature combination in the stress field abnormal area table, the potential instability mode (such as shear sliding, liquefaction) of each abnormal area is identified by using a preset rule or classification model to obtain a regional instability mode table; the characteristics of the current abnormal area are similar to the characteristics of historical cases with the same mode in the disaster precursor feature library. The similarity (cosine similarity or weighted Euclidean distance) is calculated to obtain a historical similar case matching table; for the area identified as the shear sliding mode, the safety factor F_s and the change rate dF_s / dt of the potential sliding surface are extracted to obtain a sliding surface safety factor evaluation table; for the area identified as the liquefaction mode, the pore water pressure ratio r_u and material sensitivity are extracted to evaluate the liquefaction risk and obtain a liquefaction risk evaluation table; the historical similarity, sliding surface safety factor and material sensitivity are comprehensively analyzed to obtain the liquefaction risk evaluation table. Based on the results of numerical assessment and liquefaction risk assessment, the local probability of instability in each abnormal area (point) is calculated to obtain a point instability probability table; spatial interpolation (such as co-kriging) and spatial correlation analysis are performed on the point instability probability to generate a risk spatial association map of the regional instability probability distribution; high-risk areas are identified as potential triggering areas, and the cascading effects (such as impact and pore pressure diffusion) of their instability on surrounding areas (such as downstream and adjacent areas) are simulated, and additional risks are superimposed to obtain a disaster diffusion risk map that takes the diffusion effect into account; finally, combining the disaster diffusion risk map with the time change rate of key precursor parameters (such as safety factor and stress gradient), historical data are used to predict the time window when the critical instability state is reached, and ultimately a regional instability risk map that includes spatial distribution and time prediction is formed.

[0128] The present invention is therefore intended to be illustrative and non-restrictive in all respects, with the scope of the invention being defined by the appended claims rather than the foregoing description, and all changes that come within the meaning and range of equivalents of the application documents are intended to be embraced therein.

[0129] The foregoing description is intended only to provide specific embodiments of the present invention, which will enable those skilled in the art to understand and implement the present 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 present invention. Therefore, the present invention is not intended to be limited to the embodiments shown herein, but is to be construed in the widest possible manner consistent with the principles and novel features disclosed herein.

Claims

1. A reservoir operation status prediction and analysis method based on environmental data, characterized in that: The following steps are involved: Step S1: Using a vertical resistivity probe array, a low-frequency AC signal is used to collect the original resistance response sequence of the reservoir sediment at different depths; the original resistance response sequence is subjected to polarization effect compensation processing to obtain a sediment resistivity layer spectrum; Step S2: collecting reservoir environmental data; performing time-series coupling analysis on the reservoir environmental data and the sediment resistivity layer spectrum to obtain a sediment stability influencing factor spectrum; Step S3: constructing a sediment layered structure model based on the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum; performing basic stress field calculation on the sediment layered structure model to obtain an initial stress distribution map; 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 map of the bottom mud based on the effective stress distribution map; Calculate stress gradient anomaly based on the effective stress distribution map to obtain a stress gradient anomaly area map; The stress evolution of potential sliding surface is analyzed based on the stress gradient anomaly area map and the intensity ratio distribution map to obtain the stress field map between the bottom mud layers; Step S4: Evaluate the reservoir instability probability based on the inter-layer stress field map of the sediment to obtain a regional instability risk map; Based on the regional instability risk map, a critical disaster warning is carried out and a reservoir sediment disaster warning report is obtained.

2. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is 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 reservoir sediment area coordinates 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-10 Hz with a signal strength of 5-20 mA to the probe array according to the probe position depth table. Repeat the measurement five times for each depth point with a sampling rate of 100 Hz and each measurement lasting 30 seconds. Finally, obtain the original resistance response sequence. Step S13: performing polarization effect compensation on the original resistance response sequence to obtain a set of corrected resistivity values; Step S14: performing depth interpolation processing on the corrected resistivity value set to obtain a depth distribution curve; Step S15: constructing a sediment resistivity layer spectrum according to the depth distribution curve.

3. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is characterized in that: Step S2 includes the following steps: Step S21: pre-processing the reservoir environmental data to obtain a multi-source environmental parameter time series table; Step S22: performing a lag correlation analysis on the sediment resistivity layer spectrum and the multi-source environmental parameter time series table to obtain a factor-depth response matrix; Step S23: Performing environmental factor interaction analysis on the factor-depth response matrix to obtain a factor interaction effect diagram; Step S24: Extract critical stability conditions from the factor interaction effect diagram and the factor-depth response matrix to obtain a spectrum of factors affecting sediment stability.

4. The reservoir operation status prediction and analysis method based on environmental data according to claim 3 is characterized in that: The environmental factor interaction analysis in step S2 includes: Construct a factor correlation network diagram based on the factor-depth response matrix; Perform partial correlation analysis on the factor correlation network diagram to obtain the factor direct influence matrix; Classify and identify the interaction patterns of the factor direct influence matrix to obtain the factor interaction pattern table; The time-lag interaction effect is analyzed based on the factor interaction pattern table to obtain the time-lag interaction matrix; Determine a conditional threshold table based on the time-lag interaction matrix; Perform depth difference analysis on the factor-depth response matrix according to the condition threshold table to obtain the depth interaction distribution map; Draw a factor interaction effect diagram based on the depth interaction distribution diagram, factor interaction pattern table, and conditional threshold table.

5. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is characterized in that: Constructing the sediment layer structure model in step S3 includes: The sediment physical parameters are inverted based on the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum to obtain the sediment physical property profile; According to the resistivity layer spectrum of the bottom mud, the resistivity layer peak value of the bottom mud physical property profile is detected to obtain the layer interface candidate point set; Perform characteristic clustering and stratification on the candidate point set of the layer interface to obtain a physical characteristic stratification scheme; Identify the distribution map of the interlayer transition zones of the sediment physical property profile according to the physical property stratification scheme; Identify the key weak layer map of the distribution map of the transition zone between layers; Underwater terrain fusion processing is performed according to the physical characteristics layering scheme to obtain a three-dimensional layered terrain model; Perform layered parameter space interpolation on the three-dimensional layered terrain model to obtain the layered physical parameter field; The sediment layered structure model is generated based on the layered physical parameter field and key weak layer map.

6. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is characterized in that: The foundation stress field calculation in step S3 includes: The stress-strain relationship of the sediment layered structure model is calculated to obtain the sediment mechanical response parameter set; The sediment layered structure model is discretized into a finite element grid with a horizontal grid size of 1-5 meters and a vertical grid size of 5-20 centimeters. The grid density is adjusted with depth, and the key layers are denser. Then, the corresponding mechanical response parameters are assigned to each grid cell according to the sediment mechanical response parameter set, and finally the sediment finite element grid model is obtained. Calculate the initial stress distribution diagram of the sediment finite element mesh model.

7. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is characterized in that: The pore water pressure calculation in step S3 includes: The dynamic permeability coefficient field is obtained by dynamically correcting the layered permeability coefficient of the sediment resistivity layer spectrum and the sediment stability influencing factor spectrum. Perform multi-source water pressure superposition calculation based on the dynamic permeability coefficient field to obtain a pressure component analysis table; Identify the interface pressure jump according to the pressure component analysis table and obtain the pressure jump position map; Analyze the pressure wave velocity distribution diagram based on the pressure component analysis table; The critical instability pressure is predicted based on the pressure wave velocity distribution map and the pressure jump position map to obtain the pore pressure risk map; Generate a pore water pressure distribution map based on the initial stress distribution map and the pore pressure risk map; The effective stress distribution map is calculated based on the pore water pressure distribution map and the initial stress distribution map.

8. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is characterized in that: The stress gradient anomaly calculation in step S3 includes: Compute a multi-scale gradient tensor set for effective stress distribution maps; Compute stress curl profiles based on a multi-scale gradient tensor set; Calculate the critical gradient ratio map based on the multi-scale gradient tensor set and the sediment layer structure model; Perform gradient time evolution processing on the multi-scale gradient tensor set to obtain the gradient evolution trend graph; A stress gradient anomaly area map is generated based on the critical gradient ratio map, stress curl distribution map and gradient evolution trend map.

9. The reservoir operation status prediction and analysis method based on environmental data according to claim 1 is characterized in that: The potential sliding surface stress evolution analysis in step S3 includes: Identify potential sliding surfaces based on the stress gradient anomaly area map and the intensity ratio distribution map to obtain a potential sliding surface distribution map; Calculate the sliding surface safety factor table for the potential sliding surface distribution map; The temporal and spatial evolution of the stress field is analyzed based on the effective stress distribution diagram, the sliding surface safety factor table and the stability influencing factor spectrum, and the stress field spectrum between the bottom mud layers is obtained.

10. The reservoir operation status prediction and analysis method based on environmental data according to claim 1, characterized in that: Step S4 includes the following steps: Step S41: Acquire historical disaster data of reservoir sediment; extract historical disaster features from the historical disaster data of reservoir sediment to obtain a disaster precursor feature library; Step S42: performing stress field anomaly detection based on the stress field map between sediment layers to obtain a stress field anomaly area table; Step S43: performing an instability probability assessment based on the stress field anomaly area table and the disaster precursor feature library to obtain a regional instability risk map; Step S44: Calculate the critical index based on the regional instability risk map and the stress field abnormal area table to obtain the sediment disaster critical index; Step S45: generating and issuing warning information based on the sediment disaster critical index to obtain a reservoir sediment disaster warning report.

Images