Method for predicting slope catastrophe process coupled with cumulative damage mechanism

Abstract

This invention relates to the field of computer data processing, specifically to a method for predicting slope disaster processes using coupled cumulative damage mechanisms. The method involves acquiring multi-source time-series monitoring data of the slope; mapping this data to a three-dimensional spatial feature field through radial basis function interpolation; constructing a physical information neural network by embedding the slope soil and rock damage mechanics constitutive equation and the cumulative damage evolution equation as partial differential constraint terms into the loss function of the neural network; inputting the three-dimensional spatial feature field into the network, calculating the cumulative damage tensor of each mass point within the slope space using forward modeling, and simultaneously minimizing the data prediction residual and the physical equation residual during backpropagation to update the network parameters; and determining the disaster occurrence time based on the cumulative damage tensor combined with the strain energy release rate criterion. This invention improves long-term prediction generalization ability and physical interpretability by constraining the data fitting process through physical equations, and reduces dependence on labeled data.

Description

Technical Field

[0001] This invention relates to the field of computer data processing, and more specifically to a method for predicting slope disaster processes by coupling cumulative damage mechanisms. Background Technology

[0002] Current slope disaster prediction methods largely employ purely data-driven machine learning models. These models acquire multi-source time-series monitoring data, such as slope surface displacement, deep strain, and pore water pressure, and directly use this data as input features. They then utilize deep neural networks to establish an end-to-end mapping between the monitoring data and the slope stability state. During the network training phase, the network weight parameters are iteratively updated solely by reducing the numerical error between the network's predicted output and historical labeled data. The entire fitting process treats the slope soil and rock mass as a black box system, relying solely on data distribution characteristics for function approximation, without incorporating the mechanical transmission laws and damage evolution mechanisms of the soil and rock mass.

[0003] The aforementioned data-driven prediction methods suffer from weak generalization ability and a lack of physical mechanism explanation in predicting long-term slope disaster processes. Because the fundamental mechanical process of slope instability—damage accumulating to macroscopic fracture—is not incorporated into the network training constraints, the data trajectory fitted by the network will deviate from the physically feasible region when monitoring data is outside the historical distribution or when data sparsity and noise interference exist. This leads to a sharp increase in long-term prediction errors and an inability to provide a physical quantity explanation consistent with the mechanical mechanism for the predicted disaster state. Summary of the Invention

[0004] To achieve the above objectives, the technical solution adopted by this invention is as follows: a slope disaster process prediction method coupled with a cumulative damage mechanism, comprising: acquiring multi-source time-series monitoring data of slope surface displacement, deep strain, and pore water pressure; mapping the multi-source time-series monitoring data into a three-dimensional spatial feature field through radial basis function interpolation; constructing a physical information neural network by embedding the slope rock and soil damage mechanics constitutive equation and the cumulative damage evolution equation as partial differential constraint terms into the loss function of the neural network; inputting the three-dimensional spatial feature field into the physical information neural network, and performing forward modeling to calculate the cumulative damage tensor of each mass point in the slope space, wherein the physical information neural network simultaneously minimizes the data prediction residual and the physical equation residual during backpropagation to update the network parameters; and determining the moment of macroscopic fracture surface penetration and the time of disaster occurrence based on the cumulative damage tensor output by the physical information neural network and the strain energy release rate criterion.

[0005] Preferably, mapping the multi-source time-series monitoring data into a three-dimensional spatial feature field through radial basis function interpolation includes: extracting feature vectors of different physical dimensions from the multi-source time-series monitoring data; reconstructing the initial features of the missing spatiotemporal nodes using inverse distance weighted interpolation for the missing spatiotemporal nodes in the feature vectors; constructing a radial basis interpolation matrix with Gaussian kernel function as the basis function based on the coordinates of the three-dimensional spatial grid nodes determined by slope geological survey; multiplying the reconstructed feature vectors with the radial basis interpolation matrix to project the monitoring features of the time series dimension onto the three-dimensional spatial grid nodes, eliminating the scale differences of different physical dimensions, and generating the three-dimensional spatial feature field with unified spatiotemporal resolution and integrating displacement, strain, and pore water pressure.

[0006] Preferably, a physical information neural network is constructed by embedding the slope soil and rock damage mechanics constitutive equation and the cumulative damage evolution equation as partial differential constraint terms into the loss function of the neural network. This includes: extracting the stress-strain partial derivative relationship from the slope soil and rock damage mechanics constitutive equation and the damage rate partial derivative relationship from the cumulative damage evolution equation; converting the stress-strain partial derivative relationship and the damage rate partial derivative relationship into a system of partial differential equations with respect to spatial and temporal coordinates; calculating the spatial and temporal gradients of the stress components, strain components, and damage variables predicted by the physical information neural network and substituting them into the system of partial differential equations; and adding the sum of squared residuals of the system of partial differential equations as the partial differential constraint term to the loss function.

[0007] Preferably, the physical information neural network synchronously minimizes the data prediction residual and the physical equation residual during backpropagation to update the network parameters, including: calculating the mean square error between the cumulative damage tensor of the predicted output of the physical information neural network and the actual monitoring label as the data prediction residual; calculating the sum of squared residuals of the partial differential equation system in the partial differential constraint terms as the physical equation residual; introducing dynamic weight coefficients to perform a weighted summation of the data prediction residual and the physical equation residual to construct a total loss function; and during backpropagation, performing chain-law differentiation on the partial derivatives of the weights of each layer of the physical information neural network according to the total loss function, and iteratively updating the weights of each layer using an adaptive moment estimation optimizer until the total loss function converges to a preset threshold.

[0008] Preferably, the three-dimensional spatial feature field is input into a physical information neural network to perform forward modeling calculation of the cumulative damage tensor of each mass point in the slope space. This includes: inputting the multimodal features of each grid node in the three-dimensional spatial feature field into the feature encoding layer of the physical information neural network to extract higher-order damage latent variables; concatenating the higher-order damage latent variables with the spatial and temporal coordinates of the grid nodes and inputting the result into a fully connected hidden layer to output the stress components, strain components, and damage variables of each mass point; assembling the stress components and strain components into a second-order stress tensor and a second-order strain tensor according to the Cauchy stress theorem and the definition of the strain tensor; calculating the effective damage stress in combination with the damage variables; and performing a tensor product operation between the effective damage stress and the second-order strain tensor to obtain the cumulative damage tensor of each mass point.

[0009] Preferably, based on the cumulative damage tensor output by the physical information neural network and combined with the strain energy release rate criterion, the determination of the moment of breakthrough of the macroscopic fracture surface of the slope and the time of disaster occurrence includes: extracting the main damage direction and main damage value of the cumulative damage tensor, calculating the gradient of the main damage value between adjacent spatial grids to obtain the damage localization bandwidth; performing feature decomposition on the cumulative damage tensor to obtain the damage feature vector field, integrating along the damage feature vector field to track the spatial distribution trajectory of the potential fracture surface; calculating the strain energy density integral along the spatial distribution trajectory to obtain the total strain energy release rate of the potential fracture surface; when the total strain energy release rate exceeds the fracture toughness of the soil and rock mass and the damage localization bandwidth converges to the feature thickness, recording the current time step as the moment of breakthrough of the macroscopic fracture surface of the slope and outputting the time of disaster occurrence.

[0010] Preferably, constructing a radial basis interpolation matrix with a Gaussian kernel function as the basis function includes: obtaining the distribution parameters of the elastic modulus of the soil and rock mass and the geometric parameters of the geological discontinuities in the slope geological survey; calculating the local stiffness coefficient based on the distribution parameters of the elastic modulus of the soil and rock mass, and calculating the spatial attenuation coefficient based on the geometric parameters of the geological discontinuities; using the product of the local stiffness coefficient and the spatial attenuation coefficient as the width parameter of the Gaussian kernel function, wherein the width parameter controls the spatial correlation strength between different grid nodes in the radial basis interpolation matrix; and solving the exponential term of the Gaussian kernel function based on the Euclidean distance between each grid node and the width parameter to construct a non-stationary radial basis interpolation matrix in which the width parameter varies with spatial location.

[0011] Preferably, the spatial and temporal gradients of the stress components, strain components, and damage variables predicted by the physical information neural network and substituted into the partial differential equations include: using automatic differentiation technology to calculate the partial derivatives of the stress components output by the physical information neural network with respect to spatial coordinates to obtain the stress gradient; calculating the partial derivatives of the strain components with respect to spatial coordinates to obtain the strain gradient; and calculating the partial derivatives of the damage variables with respect to time coordinates to obtain the damage evolution rate; substituting the stress gradients into the three-dimensional spatial equilibrium differential equations to construct the force equilibrium residuals; substituting the strain gradients into the geometric equations to construct the deformation compatibility residuals; and substituting the damage evolution rate into the cumulative damage evolution equations to construct the damage evolution residuals; and summing the force equilibrium residuals, the deformation compatibility residuals, and the damage evolution residuals to form the comprehensive residuals of the partial differential equations.

[0012] Preferably, a total loss function is constructed by weighting and summing the data prediction residual and the physical equation residual using dynamic weighting coefficients. This includes: calculating the partial derivative norms of the data prediction residual and the physical equation residual with respect to the weights of each layer of the physical information neural network in each iteration of backpropagation; extracting the ratio between the partial derivative norm of the data prediction residual and the partial derivative norm of the physical equation residual as a gradient ratio coefficient; constructing an exponential decay adjustment function based on the gradient ratio coefficient, and using the output value of the exponential decay adjustment function as the dynamic weighting coefficient; weighting the physical equation residual using the dynamic weighting coefficient, inversely weighting the data prediction residual, and summing the two weighted residuals to obtain the total loss function.

[0013] Preferably, the spatial distribution trajectory of the potential fracture surface is tracked by integrating the damage feature vector field, including: selecting spatial grid nodes in the three-dimensional space of the slope whose main damage value exceeds the initial damage threshold as fracture initiation seed points; starting from the fracture initiation seed points, determining a local search direction based on the damage feature vector field, and selecting the adjacent grid nodes with the largest main damage value in the neighborhood of the local search direction as the next tracking nodes; fitting the fracture initiation seed points and the sequentially selected next tracking nodes with spline curves to generate the spatial distribution trajectory of the potential fracture surface; extracting the second-order stress tensor and the second-order strain tensor from each node on the spatial distribution trajectory, calculating the strain energy density of each node, and performing curve integration along the spline curve on the strain energy density to obtain the total strain energy release rate.

[0014] Compared with existing technologies, the beneficial effects of this invention are as follows: 1. This invention constructs a physical information neural network by embedding the constitutive equation of slope soil and rock damage mechanics and the cumulative damage evolution equation as partial differential constraint terms into the loss function of the neural network. During backpropagation, the network parameters are updated by simultaneously minimizing the data prediction residual and the physical equation residual, ensuring that the network follows the solid mechanics damage evolution law while fitting the monitoring data. In scenarios with sparse data samples or distributed extrapolation, this mechanism forces the network prediction results back to the physically feasible region by constraining the physical equation residual, avoiding the divergence of the pure data model trajectory and improving the generalization ability of long-term catastrophic evolution trajectory prediction. At the same time, since the partial differential constraint terms explicitly include the stress-strain partial derivative relationship and the damage rate partial derivative relationship, the output cumulative damage tensor has interpretability at the mechanical mechanism level.

[0015] 2. This invention maps multi-source time-series monitoring data into a three-dimensional spatial feature field through radial basis function interpolation. It reconstructs missing spatiotemporal nodes and constructs a radial basis interpolation matrix using inverse distance weighted interpolation, eliminating scale differences between different physical dimensions and generating feature inputs with uniform spatiotemporal resolution, thus reducing the fitting difficulty of the network in processing multimodal data. Based on the output cumulative damage tensor, the invention determines the connection time of the macroscopic fracture surface of the slope and the time of the disaster by combining the strain energy release rate criterion. By extracting the main damage direction and main damage value, the invention obtains the damage localization bandwidth, tracks the spatial distribution trajectory of the potential fracture surface by integrating along the damage feature vector field, and calculates the total strain energy release rate. This transforms the abstract damage tensor into a specific physical criterion for fracture surface connection, realizing a logical closed loop from microscopic damage accumulation to macroscopic fracture surface evolution. Attached Figure Description

[0016] Figure 1 This is a flowchart of the multi-source time-series monitoring data acquisition and synchronous preprocessing process of the present invention.

[0017] Figure 2 This is a flowchart of the three-dimensional spatial feature field mapping of multi-source time-series monitoring data according to the present invention.

[0018] Figure 3 This is a flowchart of the physical information neural network partial differential constraint embedding construction process of the present invention.

[0019] Figure 4 This is a flowchart illustrating the dynamic update of backpropagation parameters in the physical information neural network of the present invention.

[0020] Figure 5 This is a flowchart of the forward modeling calculation of the cumulative damage tensor of slope mass points according to the present invention.

[0021] Figure 6 This is a flowchart illustrating the determination of the macroscopic fracture surface penetration and disaster occurrence time of the slope according to the present invention. Detailed Implementation

[0022] In one embodiment, the slope disaster prediction method coupled with a cumulative damage mechanism runs on a computing system comprising a data acquisition module, a data processing module, a model training module, and a disaster prediction module. The data acquisition module communicates with surface displacement monitoring stations, deep strain sensor arrays, and pore water pressure monitoring wells deployed at the slope site, acquiring multi-source time-series monitoring data in real time. The data processing module preprocesses and spatially maps the acquired raw monitoring data to generate a three-dimensional spatial feature field. The model training module constructs and trains a physical information neural network, embedding the slope soil and rock mass damage mechanics constitutive equation and the cumulative damage evolution equation as partial differential constraint terms into the network loss function. The disaster prediction module uses the trained physical information neural network to perform forward modeling calculations of the cumulative damage tensor of each mass point within the slope space, and combines this with the strain energy release rate criterion to determine the moment of macroscopic fracture surface penetration and the time of disaster occurrence.

[0023] refer to Figure 1 In one embodiment, multi-source time-series monitoring data of slope surface displacement, deep strain, and pore water pressure are acquired. Surface displacement monitoring data is acquired via a Global Navigation Satellite System (GNSS) receiver array, with each receiver recording the three-dimensional coordinate changes of monitoring points on the slope surface at fixed time intervals. Deep strain monitoring data is acquired using fiber Bragg grating strain sensors embedded at different depths on the slope. The sensors are arranged at equal intervals along the borehole axis to record the axial and lateral strains of the soil and rock at different depths. Pore water pressure monitoring data is acquired using vibrating wire pore water pressure gauges embedded below the groundwater level on the slope to record the dynamic changes in pore water pressure at different depths. All monitoring data are synchronized according to a unified timestamp, forming a multi-source time-series dataset containing time, space, and physical quantity dimensions.

[0024] refer to Figure 2 The multi-source time-series monitoring data is mapped into a three-dimensional spatial feature field through radial basis function interpolation. Feature vectors of different physical dimensions are extracted from the multi-source time-series monitoring data. These feature vectors include the x, y, and z components of surface displacement, the axial and lateral components of deep strain, and the pore water pressure scalar. For missing spatiotemporal nodes in the feature vectors, inverse distance weighted interpolation is used to reconstruct the initial features of these missing spatiotemporal nodes. Based on the coordinates of the three-dimensional spatial grid nodes determined by slope geological survey, a radial basis interpolation matrix with a Gaussian kernel function as the basis function is constructed. The reconstructed feature vectors are multiplied by the radial basis interpolation matrix, projecting the monitoring features of the time series dimension onto the three-dimensional spatial grid nodes. This eliminates scale differences between different physical dimensions, generating a three-dimensional spatial feature field with uniform spatiotemporal resolution that integrates displacement, strain, and pore water pressure.

[0025] refer to Figure 3A physical information neural network (PEN) is constructed by embedding the constitutive equations of slope soil and rock damage mechanics and the cumulative damage evolution equations as partial differential constraints into the loss function of the neural network. The stress-strain partial derivatives in the constitutive equations and the damage rate partial derivatives in the cumulative damage evolution equations are extracted. These relationships are then transformed into a system of partial differential equations with respect to spatial and temporal coordinates. The spatial and temporal gradients of the stress, strain, and damage variables predicted by the PSN are calculated and substituted into the system of partial differential equations. The sum of squared residuals from the system of partial differential equations is added as a partial differential constraint to the loss function.

[0026] refer to Figure 4 and Figure 5 The three-dimensional spatial feature field is input into the physical information neural network (PIN) to perform forward modeling of the cumulative damage tensor of each mass point within the slope space. The multimodal features of each grid node in the three-dimensional spatial feature field are input into the feature encoding layer of the PIN to extract higher-order damage latent variables. The higher-order damage latent variables are concatenated with the spatial and temporal coordinates of the grid nodes and input into a fully connected hidden layer, outputting the stress components, strain components, and damage variables of each mass point. Based on Cauchy's stress theorem and the definition of the strain tensor, the stress and strain components are assembled into a second-order stress tensor and a second-order strain tensor. The effective damage stress is calculated by combining the damage variables, and the effective damage stress is multiplied by the second-order strain tensor to obtain the cumulative damage tensor of each mass point. During backpropagation, the PIN synchronously minimizes the data prediction residual and the physical equation residual to update the network parameters. The mean square error between the cumulative damage tensor predicted by the PIN and the actual monitoring label is calculated as the data prediction residual. The sum of squared residuals of the partial differential equation system in the partial differential constraint terms is calculated as the physical equation residual. A dynamic weighting coefficient is introduced to construct the total loss function by weighted summation of the data prediction residuals and the physical equation residuals. During backpropagation, the partial derivatives of the weights of each layer in the physical information neural network are chained according to the total loss function, and the weights of each layer are iteratively updated using an adaptive moment estimation optimizer until the total loss function converges to a preset threshold.

[0027] refer to Figure 6Based on the cumulative damage tensor output by the physical information neural network, and combined with the strain energy release rate criterion, the time of breakthrough of the macroscopic fracture surface of the slope and the time of disaster occurrence are determined. The principal damage direction and principal damage value of the cumulative damage tensor are extracted, and the gradient of the principal damage value between adjacent spatial grids is calculated to obtain the damage localization bandwidth. Eigenvalue decomposition is performed on the cumulative damage tensor to obtain the damage feature vector field, and the spatial distribution trajectory of the potential fracture surface is traced by integrating along the damage feature vector field. The strain energy density integral along the spatial distribution trajectory is calculated to obtain the total strain energy release rate of the potential fracture surface. When the total strain energy release rate exceeds the fracture toughness of the soil and rock mass and the damage localization bandwidth converges to the characteristic thickness, the current time step is recorded as the time of breakthrough of the macroscopic fracture surface of the slope, and the time of disaster occurrence is output.

[0028] In this embodiment, the physical information neural network adopts a fully connected feedforward neural network structure, comprising an input layer, a feature encoding layer, multiple hidden layers, and an output layer. The input layer receives multimodal features from each grid node in the three-dimensional spatial feature field, with 6 dimensions: displacement (3D), strain (2D), and pore water pressure (1D). The feature encoding layer consists of two fully connected layers, each containing 128 neurons, using a modified linear unit (MLU) as the activation function to map the input 6-dimensional multimodal features into 128-dimensional higher-order damage latent variables. The hidden layer consists of six fully connected layers, each containing 256 neurons, using a hyperbolic tangent function as the activation function, receiving the concatenated higher-order damage latent variables, spatial coordinates, and temporal coordinates as input. The output layer contains 7 neurons, outputting the 6 independent components of the stress tensor and 1 damage variable.

[0029] In this embodiment, the slope soil and rock damage mechanics constitutive equation adopts an isotropic elastic damage constitutive model, and its expression is as follows: ;in, These are components of the second-order stress tensor. The damage variable has a range of values. , This indicates that the soil and rock mass was not damaged. This indicates that the rock and soil mass has been completely destroyed. These are the components of the fourth-order elastic stiffness tensor. These are the components of the second-order strain tensor.

[0030] For isotropic elastic materials, the fourth-order elastic stiffness tensor can be expressed as: ;in, and Let Lamé constant be . Let Kronecker function be used when hour ,when hour .

[0031] In this embodiment, the cumulative damage evolution equation adopts a damage evolution model based on strain energy density, and its expression is: ;in, For damage evolution rate, and For material constants, For strain energy release rate, The damage initiation threshold, Let be the Herveside step function, which has a value of 1 when the independent variable is greater than or equal to 0, and a value of 0 when the independent variable is less than 0.

[0032] strain energy release rate The expression is: .

[0033] In this embodiment, the expression for the three-dimensional spatial equilibrium differential equation is: ;in, For spatial coordinate components, Density of the rock and soil mass For volume force components, in a gravitational field , This is the acceleration due to gravity.

[0034] In this embodiment, the expression for the geometric equation is: ;in, These are the components of the displacement vector.

[0035] In this embodiment, the expression for the total loss function is: ;in, For the total loss function, For dynamic weighting coefficients, To predict residuals from the data, This represents the residual of the physical equation.

[0036] Data prediction residuals The expression is: ;in, The number of training samples. The first prediction of the physical information neural network The damage variable values ​​for each sample, For the first The true damage variable label value of each sample.

[0037] Physical equation residuals The expression is: ;in, For force balance residuals, To ensure deformation-compatible residuals, This refers to the damage evolution residual.

[0038] In this embodiment, the expression for the cumulative damage tensor is: ;in, These are the components of the second-order cumulative damage tensor. The equivalent stress is expressed as follows: ;in, These are components of the deviatoric stress tensor. .

[0039] In this embodiment, the expression for strain energy density is: .

[0040] The expression for the total strain energy release rate is: ;in, For the spatial distribution trajectory of the potential fracture surface, This is the unit normal vector of the fracture surface.

[0041] In this embodiment, the physical information neural network structure parameters are set as shown in Table 1.

[0042] Table 1. Parameter settings for the physical information neural network structure

[0043]

[0044] In Table 1, the input layer receives 6-dimensional multimodal features, including 3 components of displacement, 2 components of strain, and 1 component of pore water pressure. The feature encoding layer consists of two fully connected layers with ReLU activation, mapping the input features to 128-dimensional higher-order damage latent variables. The concatenation layer concatenates the 128-dimensional higher-order damage latent variables with 3-dimensional spatial coordinates and 1-dimensional time coordinates to form a 132-dimensional input vector. All six hidden layers use the Tanh activation function, each containing 256 neurons, to learn complex nonlinear mechanical relationships. The output layer outputs 7 physical quantities, including 6 independent components of the stress tensor and 1 damage variable.

[0045] In this embodiment, the training dataset consists of historical monitoring data and numerical simulation data. The historical monitoring data comes from the past five years of on-site slope monitoring records, including surface displacement, deep strain, and pore water pressure data, as well as corresponding slope stability state labels. The numerical simulation data is generated using finite element software to simulate the damage evolution process of the slope under different working conditions, generating a large amount of stress, strain, and damage data. The training dataset is divided into training, validation, and test sets in a 7:2:1 ratio. During training, the learning rate of the adaptive moment estimation optimizer is set to 0.001, the weight decay coefficient is set to 0.0001, the batch size is set to 32, and the maximum number of iterations is set to 10,000. Training is terminated early when the validation set loss no longer decreases after 500 consecutive iterations.

[0046] In this embodiment, the three-dimensional spatial mesh of the slope is divided using tetrahedral elements, with 10,000 mesh nodes and 50,000 elements. The mesh resolution is adjusted according to the slope geological conditions and the distribution of monitoring points. Mesh refinement is performed in areas with dense monitoring points and near geological discontinuities to improve computational accuracy in these areas. The temporal resolution of the three-dimensional spatial feature field is consistent with the sampling interval of the monitoring data, which is 1 hour.

[0047] This embodiment embeds the constitutive equation of slope soil and rock damage mechanics and the cumulative damage evolution equation as partial differential constraints into the loss function of a neural network. During backpropagation, it simultaneously minimizes the data prediction residual and the physical equation residual to update the network parameters, ensuring that the network follows the solid mechanics damage evolution law while fitting the monitoring data. Radial basis function interpolation maps multi-source time-series monitoring data into a three-dimensional spatial feature field, eliminating scale differences between different physical dimensions and generating feature inputs with uniform spatiotemporal resolution. Based on the output cumulative damage tensor and combined with the strain energy release rate criterion, the timing of the macroscopic fracture surface penetration and the time of catastrophic event are determined, realizing a logical closed loop from microscopic damage accumulation to macroscopic fracture surface evolution.

[0048] In a preferred embodiment, feature vectors of different physical dimensions are extracted from the multi-source time-series monitoring data. For missing spatiotemporal nodes in the feature vectors, the initial features of the missing spatiotemporal nodes are reconstructed using inverse distance weighted interpolation. During the acquisition process, multi-source time-series monitoring data may experience data loss due to sensor malfunctions, communication interruptions, or environmental interference, resulting in incomplete spatiotemporal feature vectors. For each missing spatiotemporal node, valid monitoring data points within a certain spatial and temporal range are searched, and the initial feature value of the missing node is calculated using the inverse distance weighted interpolation method. The basic idea of ​​inverse distance weighted interpolation is that points closer to each other have a greater influence on the interpolation point, and the weight is inversely proportional to the power of the distance.

[0049] Based on the coordinates of the three-dimensional spatial grid nodes determined by the slope geological survey, a radial basis function (RBF) interpolation matrix is ​​constructed using the Gaussian kernel function as the basis function. The elastic modulus distribution parameters of the soil and rock mass and the geometric parameters of the geological discontinuities in the slope geological survey are obtained. Local stiffness coefficients are calculated based on the elastic modulus distribution parameters of the soil and rock mass, and spatial attenuation coefficients are calculated based on the geometric parameters of the geological discontinuities. The product of the local stiffness coefficient and the spatial attenuation coefficient is used as the width parameter of the Gaussian kernel function, where the width parameter controls the spatial correlation strength between different grid nodes in the RBF interpolation matrix. Based on the Euclidean distance between each grid node and the width parameter, the exponential term of the Gaussian kernel function is solved, and a non-stationary RBF interpolation matrix with a width parameter varying with spatial location is constructed.

[0050] The reconstructed feature vectors are multiplied by the radial basis interpolation matrix, projecting the time-series monitoring features onto the 3D spatial grid nodes. This eliminates scale differences between different physical dimensions, generating a 3D spatial feature field with uniform spatiotemporal resolution that integrates displacement, strain, and pore water pressure. For each time step, the reconstructed feature vectors are multiplied by the radial basis interpolation matrix to obtain the feature values ​​of all 3D spatial grid nodes at that time step. Feature values ​​with different physical dimensions are standardized to map them to the same numerical range, eliminating the impact of dimensional differences on subsequent neural network training. The 3D spatial grid feature values ​​from all time steps are combined to form a 3D spatial feature field with uniform spatiotemporal resolution.

[0051] In this embodiment, the expression for inverse distance weighted interpolation is: ;in, Missing node The interpolation results, The number of valid data points participating in the interpolation. For the first Feature values ​​of each valid data point For the first The weights of each valid data point.

[0052] Weight The expression is: ;in, Missing node With the 1 valid data point The spatiotemporal distance between them It is a power exponent.

[0053] Spacetime Distance The expression is: ;in, and These are the spatial coordinates of the missing node and the valid data point, respectively. and These are the timestamps for missing nodes and valid data points, respectively. This is the conversion factor between time distance and spatial distance.

[0054] In this embodiment, the expression for the Gaussian kernel function is: ;in, The Euclidean distance between two grid nodes. is the width parameter of the Gaussian kernel function.

[0055] Radial basis interpolation matrix elements for: ;in, For the first The grid node and the first Euclidean distance between grid nodes For the first Width parameters at each grid node.

[0056] In this embodiment, the local stiffness coefficient The expression is: ;in, For the first Elastic modulus of soil and rock at each grid node This represents the average elastic modulus of the slope's soil and rock mass.

[0057] Spatial attenuation coefficient The expression is: ;in, For the first The distance from each grid node to the nearest geological discontinuity. This refers to the influence range parameter of the geological discontinuity surface.

[0058] Width parameter of Gaussian kernel function The expression is: ;in, This is the baseline width parameter.

[0059] In this embodiment, the standardization of different physical dimensions is performed using the Z-score standardization method, the expression of which is: ;in, These are the standardized eigenvalues. These are the original eigenvalues. This is the mean of the physical quantity across all spatiotemporal nodes. This represents the standard deviation of the physical quantity across all spatiotemporal nodes.

[0060] In this embodiment, the standardized parameters of different physical quantities are shown in Table 2.

[0061] Table 2 Standardized parameters for different physical quantities

[0062]

[0063] In Table 2, these parameters were obtained through statistical analysis of historical monitoring data. Z-score standardization maps the feature values ​​of different physical dimensions to a standard normal distribution with a mean of 0 and a standard deviation of 1, eliminating the impact of dimensional differences on neural network training and improving the network's convergence speed and prediction accuracy.

[0064] In this embodiment, the power exponent of inverse distance weighted interpolation Set to 2, the conversion factor between time distance and spatial distance. The time frame is set to 0.1 hours / meter. The spatial search radius for valid data points used in the interpolation is set to 50 meters, and the time search window is set to 24 hours. When the number of valid data points for a missing node within the search range is less than 3, the historical average value of that physical quantity at that spatial location is used as the interpolation result.

[0065] In this embodiment, geological discontinuities include faults, joints, fissures, and weak interlayers. The geometric parameters of these discontinuities are obtained through slope geological surveys, including their strike, dip, dip angle, length, width, and spacing. The parameters defining the influence range of the geological discontinuities are also included. The determination is based on the size and nature of the discontinuity; for large faults, Set to 20 meters; for medium-sized joints, Set to 10 meters; for small cracks, Set to 5 meters.

[0066] In this embodiment, the coordinates of the three-dimensional spatial grid nodes are generated using a slope digital elevation model and a geological profile. First, surface grid nodes are generated based on the slope digital elevation model. Then, internal grid nodes are generated based on the geological profile. Finally, a complete three-dimensional spatial grid is generated by dividing the grid using tetrahedral elements. The grid resolution is set to 1 meter near the slope surface and potential sliding surface, 5 meters inside the slope, and 10 meters deep within the slope.

[0067] This embodiment reconstructs the initial features of missing spatiotemporal nodes using inverse distance weighted interpolation, solving the problem of missing monitoring data and ensuring data integrity. By constructing a non-stationary radial basis interpolation matrix with a width parameter varying with spatial location, the influence of the elastic modulus distribution of soil and rock and the spatial correlation strength of geological discontinuities is considered, improving the accuracy of spatial interpolation. Standardization eliminates scale differences between different physical dimensions, generating a three-dimensional spatial feature field with uniform spatiotemporal resolution, providing high-quality input data for the subsequent training of the physical information neural network.

[0068] In a preferred embodiment, the stress-strain partial derivative relationship in the constitutive equation of slope soil and rock damage mechanics and the damage rate partial derivative relationship in the cumulative damage evolution equation are extracted. The constitutive equation of slope soil and rock damage mechanics describes the relationship between stress and strain, including the partial derivative of stress with respect to strain and the partial derivative of stress with respect to the damage variable. The cumulative damage evolution equation describes the change of the damage variable over time, including the partial derivative of the damage variable with respect to time and the partial derivative of the damage variable with respect to the strain energy release rate.

[0069] The stress-strain partial derivative relationship and the damage rate partial derivative relationship are transformed into a system of partial differential equations with respect to spatial and temporal coordinates. The stress-strain partial derivative relationship, combined with the three-dimensional equilibrium differential equations and geometric equations, forms a system of partial differential equations concerning displacement, stress, and strain. The damage rate partial derivative relationship, combined with the expression for the strain energy release rate, forms a system of partial differential equations concerning damage variables. These partial differential equations are combined to form a complete system of partial differential equations describing the damage evolution process of slope soil and rock.

[0070] The stress gradient is obtained by calculating the partial derivatives of the stress component with respect to spatial coordinates from the output of the physical information neural network using automatic differentiation technology. The strain gradient is obtained by calculating the partial derivatives of the strain component with respect to spatial coordinates, and the damage evolution rate is obtained by calculating the partial derivatives of the damage variable with respect to time coordinates. Automatic differentiation technology automatically calculates the partial derivatives of the function with respect to the input variables using the chain rule, enabling efficient and accurate calculation of the partial derivatives of the neural network output with respect to spatial and time coordinates.

[0071] The stress gradient is substituted into the three-dimensional equilibrium differential equation to construct the force equilibrium residual; the strain gradient is substituted into the geometric equation to construct the deformation compatibility residual; and the damage evolution rate is substituted into the cumulative damage evolution equation to construct the damage evolution residual. The force equilibrium residual, deformation compatibility residual, and damage evolution residual are summed to form the comprehensive residual of the partial differential equation system. The sum of squared residuals of the partial differential equation system is added as a partial differential constraint term to the loss function.

[0072] The mean squared error between the cumulative damage tensor predicted by the physical information neural network and the actual monitoring label is calculated as the data prediction residual. The sum of squared residuals of the partial differential equation system in the partial differential constraint terms is calculated as the physical equation residual. In each iteration of backpropagation, the partial derivative norms of the data prediction residual and the physical equation residual with respect to the weights of each layer of the physical information neural network are calculated. The ratio between the partial derivative norm of the data prediction residual and the partial derivative norm of the physical equation residual is extracted as the gradient ratio coefficient. An exponential decay adjustment function is constructed based on the gradient ratio coefficient, and the output value of the exponential decay adjustment function is used as the dynamic weight coefficient. The physical equation residual is weighted by the dynamic weight coefficient, and the data prediction residual is inversely weighted. The sum of the two weighted residuals is used to obtain the total loss function.

[0073] During the backpropagation process, the partial derivatives of the weights of each layer of the physical information neural network are chained according to the total loss function, and the weights of each layer are iteratively updated using an adaptive moment estimation optimizer until the total loss function converges to a preset threshold.

[0074] The multimodal features of each grid node in the three-dimensional spatial feature field are input into the feature encoding layer of the physical information neural network to extract higher-order damage latent variables. The higher-order damage latent variables are concatenated with the spatial and temporal coordinates of the grid nodes and then input into a fully connected hidden layer to output the stress components, strain components, and damage variables of each mass point. Based on Cauchy's stress theorem and the definition of the strain tensor, the stress and strain components are assembled into a second-order stress tensor and a second-order strain tensor. The effective damage stress is calculated using the damage variables, and the effective damage stress is multiplied by the second-order strain tensor to obtain the cumulative damage tensor of each mass point.

[0075] The principal damage direction and principal damage value of the cumulative damage tensor are extracted, and the gradient of the principal damage value between adjacent spatial grids is calculated to obtain the damage localization bandwidth. Eigenvalue decomposition is performed on the cumulative damage tensor to obtain a damage feature vector field, and the spatial trajectory of the potential fracture surface is tracked by integrating along the damage feature vector field. Spatial grid nodes in the three-dimensional space of the slope where the principal damage value exceeds the initial damage threshold are selected as fracture initiation seed points. Starting from the fracture initiation seed points, a local search direction is determined based on the damage feature vector field, and the adjacent grid node with the largest principal damage value in the neighborhood of the local search direction is selected as the next tracking node. Spline curve fitting is performed between the fracture initiation seed points and the sequentially selected next tracking nodes to generate the spatial trajectory of the potential fracture surface.

[0076] The second-order stress tensor and the second-order strain tensor are extracted from each node on the spatial distribution trajectory. The strain energy density of each node is calculated, and the total strain energy release rate is obtained by curve integration along the spline curve. When the total strain energy release rate exceeds the fracture toughness of the soil and rock mass and the damage localization bandwidth converges to the characteristic thickness, the current time step is recorded as the moment when the macroscopic fracture surface of the slope is penetrated, and the time of the catastrophe is output.

[0077] In this embodiment, the automatic differentiation technique employs forward-mode automatic differentiation. For each physical quantity output by the neural network, a computational graph of it with respect to spatial and temporal coordinates is constructed, and then partial derivatives are calculated sequentially from input to output using the chain rule. Stress gradient strain gradient and damage evolution rate All were calculated using automatic differentiation techniques.

[0078] In this embodiment, the force balance residual The expression is: ;in, The number of nodes in the three-dimensional spatial mesh. For the first Stress tensor components at each grid node.

[0079] Deformation Coordination Residual The expression is: ;in, For the first Strain tensor components at each grid node For the first Displacement vector components at each grid node.

[0080] Damage evolution residual The expression is: ;in, For the first Damage variables at each grid node For the first Strain energy release rate at each grid node.

[0081] In this embodiment, the norm of the partial derivative of the data prediction residual The norm of the partial derivatives of the physical equation residuals They are respectively.

[0082] .

[0083] .

[0084] in, The set of all parameters of the neural network, including weights. and bias , The number of layers in the neural network. For the first The number of neurons in a layer.

[0085] Gradient ratio coefficient The expression is: .

[0086] The expression for the exponential decay adjustment function is: ;in, For adjustment coefficients, The target gradient ratio.

[0087] In this embodiment, the eigenvalue decomposition expression of the cumulative damage tensor is: ;in, For the first Each primary damage value satisfies , For the corresponding to the first The unit feature vector of each main damage value, i.e., the main damage direction.

[0088] Damage localization bandwidth The expression is: ;in, The gradient magnitude of the main damage value.

[0089] In this embodiment, cubic B-spline curves are used for spline curve fitting. For a set of ordered control points... The expression for a cubic B-spline curve is: ;in, For curve parameters, It is a cubic B-spline basis function.

[0090] In this embodiment, the values ​​of the mechanical parameters of the soil and rock mass are shown in Table 3.

[0091] Table 3 Values ​​of Mechanical Parameters of Soil and Rock Mass

[0092]

[0093] In Table 3, these parameters were obtained through indoor rock mechanics tests and field geological surveys. The elastic modulus and Poisson's ratio were determined through uniaxial compression tests, density through specific gravity tests, and the Lamé constant was calculated from the elastic modulus and Poisson's ratio. Damage initiation threshold and material constants are also included. and The fracture toughness was determined by triaxial compression and damage evolution tests. The characteristic thickness was determined by three-point bending tests. The characteristic thickness was determined based on the particle size distribution and structural characteristics of the slope soil and rock mass.

[0094] In this embodiment, the parameters of the adaptive moment estimation optimizer are set as follows: initial learning rate of 0.001, learning rate decay coefficient of 0.9, exponential decay rate of first-order moment estimation of 0.9, exponential decay rate of second-order moment estimation of 0.999, and numerical stability constant of . The convergence preset threshold for the total loss function is set to... When the change in the total loss function over 1000 consecutive iterations is less than At that time, it was assumed that the network training had converged.

[0095] In this embodiment, the initial damage threshold is set to 0.1. When the principal damage value of a grid node exceeds 0.1, it is marked as a seed point for fracture initiation. The local search direction is the direction of maximum principal damage, and the neighborhood range is a 3×3×3 grid cube. In each search step, the neighboring grid node with the largest principal damage value within the neighborhood of the local search direction is selected as the next tracking node. Tracking stops when the principal damage values ​​of three consecutive tracking nodes are less than 0.05.

[0096] In this embodiment, the method for predicting the time of disaster occurrence is as follows: after the macroscopic rupture surface of the slope is detected to be connected, the time required for the rupture surface to connect and for the slope to become unstable as a whole is calculated based on the expansion rate of the rupture surface and the geometric dimensions of the slope. This time is added to the moment when the macroscopic rupture surface connects to obtain the time of disaster occurrence. The expansion rate of the rupture surface is determined based on historical monitoring data and numerical simulation results, and the value ranges from 0.1 to 10 meters per hour.

[0097] This embodiment efficiently and accurately calculates the residuals of the physical equations using automatic differentiation technology. The constitutive equations of slope soil and rock damage mechanics and the cumulative damage evolution equations are embedded as partial differential constraints in the loss function of the neural network, ensuring that the network's prediction results follow the laws of solid mechanics. By introducing dynamic weight coefficients to adaptively adjust the weights of the data prediction residuals and the physical equation residuals in the total loss function, the importance of data fitting and physical constraints is balanced, improving the network's generalization ability. By tracing the spatial distribution trajectory of potential fracture surfaces along the damage feature vector field integral and combining this with the strain energy release rate criterion to determine the moment of macroscopic fracture surface penetration and the time of catastrophic event, a complete prediction process from microscopic damage accumulation to macroscopic catastrophic event occurrence is achieved.

Claims

1. A method for predicting slope catastrophic processes using coupled cumulative damage mechanisms, characterized in that, include: Acquire multi-source time-series monitoring data of slope surface displacement, deep strain, and pore water pressure; The multi-source time-series monitoring data is mapped into a three-dimensional spatial feature field through radial basis function interpolation; The physical information neural network is constructed by embedding the constitutive equation of slope rock and soil damage mechanics and the cumulative damage evolution equation as partial differential constraint terms into the loss function of the neural network. The three-dimensional spatial feature field is input into the physical information neural network to calculate the cumulative damage tensor of each mass point in the slope space in forward modeling. The physical information neural network simultaneously minimizes the data prediction residual and the physical equation residual during back propagation to update the network parameters. Based on the cumulative damage tensor output by the physical information neural network, and combined with the strain energy release rate criterion, the time of breakthrough of the macroscopic fracture surface of the slope and the time of disaster occurrence are determined.

2. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 1, characterized in that, Mapping the multi-source time-series monitoring data into a three-dimensional spatial feature field through radial basis function interpolation includes: extracting feature vectors of different physical dimensions from the multi-source time-series monitoring data; and reconstructing the initial features of the missing spatiotemporal nodes using inverse distance weighted interpolation for the missing spatiotemporal nodes in the feature vectors. Based on the coordinates of the three-dimensional spatial grid nodes determined by the slope geological survey, a radial basis interpolation matrix with the Gaussian kernel function as the basis function is constructed. The reconstructed feature vector is multiplied by the radial basis interpolation matrix, and the monitoring features of the time series dimension are projected onto the three-dimensional spatial grid nodes to eliminate the scale differences of different physical dimensions and generate the three-dimensional spatial feature field with unified spatiotemporal resolution and fusion of displacement, strain and pore water pressure.

3. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 1, characterized in that, The physical information neural network is constructed by embedding the constitutive equation of slope rock and soil damage mechanics and the cumulative damage evolution equation as partial differential constraint terms into the loss function of the neural network. This includes: extracting the stress-strain partial derivative relationship in the constitutive equation of slope rock and soil damage mechanics and the damage rate partial derivative relationship in the cumulative damage evolution equation. The stress-strain partial derivative relationship and the damage rate partial derivative relationship are converted into a system of partial differential equations with respect to spatial coordinates and time coordinates; The spatial and temporal gradients of the stress components, strain components, and damage variables predicted by the physical information neural network are calculated and substituted into the partial differential equation system. The sum of squared residuals of the partial differential equation system is then added to the loss function as the partial differential constraint term.

4. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 3, characterized in that, The physical information neural network synchronously minimizes the data prediction residual and the physical equation residual during backpropagation to update the network parameters, including: calculating the mean square error between the cumulative damage tensor of the predicted output of the physical information neural network and the actual monitoring label as the data prediction residual; The sum of squared residuals of the partial differential equation system in the partial differential constraint terms is calculated as the residual of the physical equation; A total loss function is constructed by weighting the data prediction residuals and the physical equation residuals with dynamic weighting coefficients. During backpropagation, the partial derivatives of the weights of each layer of the physical information neural network are chained according to the total loss function, and the weights of each layer are iteratively updated using an adaptive moment estimation optimizer until the total loss function converges to a preset threshold.

5. The slope disaster process prediction method based on coupled cumulative damage mechanism according to claim 2, characterized in that, The three-dimensional spatial feature field is input into the physical information neural network to perform forward modeling calculation of the cumulative damage tensor of each mass point in the slope space, including: inputting the multimodal features of each grid node in the three-dimensional spatial feature field into the feature encoding layer of the physical information neural network to extract higher-order damage latent variables; The higher-order damage hidden variables are concatenated with the spatial and temporal coordinates of the mesh nodes and then input into the fully connected hidden layer to output the stress components, strain components and damage variables of each mass point. Based on Cauchy's stress theorem and the definition of strain tensor, the stress components and strain components are assembled into a second-order stress tensor and a second-order strain tensor. The effective damage stress is calculated in combination with the damage variables. The effective damage stress and the second-order strain tensor are multiplied to obtain the cumulative damage tensor of each mass point.

6. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 5, characterized in that, Based on the cumulative damage tensor output by the physical information neural network, and combined with the strain energy release rate criterion, the time of breakthrough of the macroscopic rupture surface of the slope and the time of disaster occurrence are determined, including: extracting the main damage direction and main damage value of the cumulative damage tensor, and calculating the gradient of the main damage value between adjacent spatial grids to obtain the damage localization bandwidth. The cumulative damage tensor is decomposed to obtain a damage feature vector field, and the spatial distribution trajectory of the potential fracture surface is traced by integrating along the damage feature vector field. Calculate the strain energy density integral along the spatial distribution trajectory to obtain the total strain energy release rate of the potential fracture surface; When the total strain energy release rate exceeds the fracture toughness of the soil and rock mass and the damage localization bandwidth converges to the characteristic thickness, the current time step is recorded as the moment when the macroscopic fracture surface of the slope is connected and the time of the disaster is output.

7. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 2, characterized in that, Construct a radial basis interpolation matrix with Gaussian kernel function as the basis function, including: obtaining the distribution parameters of elastic modulus of soil and rock mass and the geometric parameters of geological discontinuities in slope geological survey; The local stiffness coefficient is calculated based on the elastic modulus distribution parameters of the soil and rock mass, and the spatial attenuation coefficient is calculated based on the geometric parameters of the geological discontinuity surface. The product of the local stiffness coefficient and the spatial attenuation coefficient is used as the width parameter of the Gaussian kernel function, wherein the width parameter controls the spatial correlation strength between different grid nodes in the radial basis interpolation matrix; Based on the Euclidean distance between each grid node and the width parameter, the exponential term of the Gaussian kernel function is solved to construct a non-stationary radial basis interpolation matrix in which the width parameter varies with spatial location.

8. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 3, characterized in that, The spatial and temporal gradients of the stress components, strain components, and damage variables predicted by the physical information neural network are obtained and substituted into the partial differential equations. This includes: using automatic differentiation technology to calculate the partial derivatives of the stress components output by the physical information neural network with respect to spatial coordinates to obtain the stress gradient; calculating the partial derivatives of the strain components with respect to spatial coordinates to obtain the strain gradient; and calculating the partial derivatives of the damage variables with respect to time coordinates to obtain the damage evolution rate. The stress gradient is substituted into the three-dimensional spatial equilibrium differential equation to construct the force equilibrium residual; the strain gradient is substituted into the geometric equation to construct the deformation compatibility residual; and the damage evolution rate is substituted into the cumulative damage evolution equation to construct the damage evolution residual. The summation of the force balance residual, the deformation coordination residual, and the damage evolution residual constitutes the comprehensive residual of the partial differential equation system.

9. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 4, characterized in that, The total loss function is constructed by weighting and summing the data prediction residual and the physical equation residual using dynamic weight coefficients, including: calculating the partial derivative norms of the data prediction residual and the physical equation residual with respect to the weights of each layer of the physical information neural network in each iteration cycle of backpropagation; The ratio between the partial derivative norm of the data prediction residual and the partial derivative norm of the physical equation residual is extracted as the gradient ratio coefficient. An exponential decay adjustment function is constructed based on the gradient ratio coefficient, and the output value of the exponential decay adjustment function is used as the dynamic weight coefficient. The physical equation residuals are weighted by the dynamic weighting coefficients, and the data prediction residuals are inversely weighted. The total loss function is obtained by summing the two weighted residuals.

10. The slope catastrophic process prediction method based on coupled cumulative damage mechanism according to claim 6, characterized in that, Tracing the spatial distribution trajectory of potential fracture surfaces by integrating the damage feature vector field includes: selecting spatial grid nodes in the three-dimensional space of the slope whose main damage value exceeds the initial damage threshold as seed points for fracture initiation; Starting from the rupture initiation seed point, a local search direction is determined based on the damage feature vector field, and the neighboring grid node with the largest main damage value in the neighborhood of the local search direction is selected as the next tracking node. The spatial distribution trajectory of the potential fracture surface is generated by fitting the fracture initiation seed point with the next selected tracking node using a spline curve. The second-order stress tensor and second-order strain tensor are extracted from each node on the spatial distribution trajectory, the strain energy density of each node is calculated, and the total strain energy release rate is obtained by curve integration of the strain energy density along the spline curve.

Images