A search method for determining the double critical slip surface of a heterogeneous soil layer slope

Abstract

The present invention relates to a search method for determining the double critical slip surface of a heterogeneous soil slope. This method uses Monte Carlo random simulation combined with genetic algorithm to generate and screen parameter samples, and determines the optimal parameter group through a multi-index similarity function and assigns it to the model. Based on the improved NSGA-II algorithm, two sets of candidate slip surfaces are generated synchronously, and a collaborative safety factor is constructed by integrating the strength reduction method and the energy method to comprehensively consider the mechanical response of the slope. The failure probability is solved through Monte Carlo simulation to judge the double critical slip surface. Finally, the optimal critical slip surface is locked by using the collaborative framework of particle swarm optimization and tabu search. This algorithm can accurately capture the double critical slip surface, is widely applicable to various geological and load conditions, improves the calculation accuracy and efficiency, and provides key data support for slope engineering design, safety assessment, disaster warning, etc.

Description

Technical Field

[0001] The present invention relates to the technical field of slope stability analysis, and particularly relates to a search method for determining a double critical slip surface of a heterogeneous soil layer slope. Background Art

[0002] Slope stability analysis is one of the core topics in the field of geotechnical engineering. The development of its theoretical system and calculation methods directly affects the reliability of engineering safety assessment and the foresight of disaster prevention and control. Traditional analysis methods are mostly based on the rigid body limit equilibrium theory framework, and the safety factor of the potential slip surface is solved by constructing static equilibrium equations. These methods simplify the sliding body into a rigid body motion and use the slice method to discretize the slip surface. The core assumption is the singularity of the slip surface morphology and the idealization of the sliding mechanism. However, actual engineering slopes are often controlled by the coupling action of multiple factors: the formation structure shows significant heterogeneity, and the weak interlayers and fracture networks form complex stress transfer paths; the spatio-temporal variability of the groundwater seepage field and geotechnical mechanical parameters leads to the strength weakening effect; the dynamic application and unloading process of external loads induce progressive failure. Under such complex working conditions, the slope instability mode often shows a composite failure mechanism of multi-slip surface cooperative action, and the single slip surface assumption is difficult to accurately characterize its mechanical essence.

[0003] Existing methods have inherent limitations at the theoretical level:

[0004] First, the ability to model the interaction effect of the slip surface is insufficient, and the energy transfer and kinematic correlation between the double slip surfaces cannot be quantified;

[0005] Second, the parameter sensitivity analysis is mostly limited to single-factor perturbation, lacking a global assessment of the system response under the co-variation of multiple parameters;

[0006] Third, the safety factor calculation criterion relies too much on the strength reduction principle and fails to integrate multiple criteria such as energy dissipation and displacement field evolution, resulting in deviations in the determination of the instability critical state.

[0007] With the development of computational mechanics and optimization algorithms, the slope stability analysis method based on numerical simulation has gradually become a research hotspot. Although numerical methods such as the finite element method and the discrete element method can more precisely depict the slope stress and strain field, they still face severe challenges in the critical slip surface search process.

[0008] Existing search methods generally have the contradiction between convergence efficiency and global optimization ability:

[0009] Deterministic algorithms (such as the gradient descent method) rely on the initial parameter setting and are prone to falling into local optimal solutions and missing the true critical slip surface;

[0010] Randomized algorithms such as genetic algorithms and simulated annealing have the potential for global exploration, but the computational cost grows exponentially when the dimensionality of the parameter space explodes. In addition, for multi-sliding surface search, a high-dimensional optimization model needs to be constructed, and it is difficult for traditional single-objective optimization frameworks to balance the dual constraints of minimizing the safety factor and maximizing the sliding volume, resulting in insufficient completeness of the solution set. At the level of grid discretization, the static uniform partitioning strategy cannot adaptively capture the strain concentration regions near the sliding surface, causing waste of computational resources and loss of accuracy.

[0011] Although existing studies have tried to introduce dynamic grid technology, the subdivision triggering mechanisms are mostly based on empirical thresholds and lack self-consistent correlations with mechanical response characteristics, making it difficult to achieve coordinated optimization of computational efficiency and accuracy. These bottlenecks restrict the interpretability of complex slope instability mechanisms, and there is an urgent need to construct a new analysis framework that integrates multi-source data-driven, intelligent optimization strategies, and physical mechanism constraints to break through the theoretical limitations and practical obstacles of traditional methods.

[0012] Traditional slope stability analysis methods face fundamental theoretical bottlenecks when dealing with the coupled instability mechanism of multi-sliding surfaces:

[0013] Firstly, the single critical sliding surface search framework simplifies the potential failure path into a single mechanical channel, ignoring the dominant role of the co-evolution of double sliding surfaces in system stability in actual engineering. Such simplifying assumptions lead to the model's inability to capture the stress redistribution, energy transfer, and kinematic correlations between sliding surfaces, resulting in systematic biases in safety factor assessment. Especially in heterogeneous slopes with weak interlayers and fault intersections, the risk of misjudgment increases significantly.

[0014] Secondly, although existing numerical optimization algorithms show certain effectiveness in single sliding surface search, their single-objective optimization mechanism is difficult to coordinate the multi-constraint contradictions among the geometric shape of the sliding surface, minimizing the safety factor, and maximizing the sliding volume, resulting in incomplete coverage of the solution set space and omission of the true critical sliding surface combinations. In addition, the static discretization strategy and empirical parameter sensitivity analysis further weaken the model's adaptability to complex geological conditions and dynamic disturbances, making it difficult to meet the requirements of high-precision engineering decisions. Summary of the Invention

[0015] The present invention provides a method for searching for double critical sliding surfaces of heterogeneous soil layer slopes, aiming to efficiently and accurately find the double critical sliding surfaces under complex working conditions of slopes and provide key data support for slope engineering design and safety assessment.

[0016] To solve the above technical problems, the present invention provides the following technical solutions:

[0017] In the first aspect, the present invention provides a method for searching for double critical sliding surfaces of heterogeneous soil layer slopes, which specifically includes the following steps:

[0018] Step 1: Obtain the geological parameters and topographic data of the area to be analyzed;

[0019] Step 2: Construct a three-dimensional real-scene model and generate a finite element numerical model according to the geological parameters and topographic data;

[0020] Step 3: Set the parameter value range for the obtained geological parameters, determine the probability distribution function of the parameters, generate and screen parameter samples using Monte Carlo stochastic simulation and genetic algorithms, and determine the optimal parameter group through a multi-index similarity function and assign it to the finite element numerical model;

[0021] Step 4: For the assigned finite element numerical model, adopt an adaptive mesh refinement technique to improve the mesh resolution of the finite element numerical model and obtain a numerical model with refined meshes;

[0022] Step 5: Based on the numerical model with refined meshes, synchronously generate two groups of sliding surface candidate sets using the improved NSGA-II algorithm, construct a collaborative safety factor by integrating the strength reduction method and the energy method, and search for parallel double sliding surfaces in the two groups of sliding surface candidate sets;

[0023] Step 6: According to the obtained parallel double sliding surfaces, solve the failure probability through Monte Carlo simulation and determine the double critical sliding surfaces;

[0024] Step 7: Based on the double critical sliding surfaces, lock the optimal critical sliding surface using the collaborative framework of particle swarm optimization and tabu search.

[0025] Furthermore, obtaining the geological parameters and topographic data includes:

[0026] Investigate the geological parameters of the geology through on-site geological exploration, including rock and soil stratification parameters, spatial distribution of soft interlayers, and dynamic data of the groundwater level;

[0027] The rock and soil stratification parameters include cohesion c, internal friction angle φ, and unit weight γ;

[0028] The spatial distribution of the soft interlayers includes thickness and dip angle;

[0029] The dynamic data of the groundwater level includes burial depth and permeability coefficient k.

[0030] Adopt UAV photogrammetry and airborne LiDAR to obtain topographic data, including digital orthophoto maps, digital elevation models DEM, and digital surface models DSM of the study area.

[0031] Furthermore, constructing a three-dimensional real-scene model and generating a finite element numerical model according to the geological parameters and topographic data includes:

[0032] Construct a 3D real - scene model from the obtained topographic data;

[0033] Import the constructed 3D real - scene model into finite - element software to generate a finite - element numerical model,

[0034] Perform an initial mesh division on the overall area in the finite - element numerical model, and discretize the entire finite - element numerical model into a finite number of elements and nodes with equal - sized grids.

[0035] Furthermore, for the obtained geological parameters, set their parameter value ranges, determine the probability distribution function of the parameters, generate and screen parameter samples using Monte Carlo random simulation and genetic algorithms, and determine the optimal parameter group through a multi - index similarity function and assign it to the finite - element numerical model, including:

[0036] Set the parameter value range and probability distribution function for the obtained geological parameters, generate parameter groups using Monte Carlo random simulation to form an initial parameter sample set. Each parameter group represents a possible parameter combination, and at the same time, embed the genetic - algorithm gene encoding, mapping the parameter group to a chromosome individual,

[0037] Randomly select M parameter groups from the initial parameter sample set to ensure that the samples cover the entire domain of the finite - element numerical model,

[0038] Evaluate the geological matching degree of each selected parameter combination through a multi - index similarity function, and screen out the optimal parameter group;

[0039] Assign values to different rock and soil layers in the finite - element numerical model according to the optimal parameter group.

[0040] Furthermore, the multi - index similarity function is:

[0041]

[0042] Among them, , represents the normalized cohesion difference, where represents the sample cohesion, represents the in - situ measured cohesion, represents the standard deviation of the in - situ cohesion, is the deviation angle of the internal friction angle, KS(GSD) is the Kolmogorov - Smirnov test statistic of the particle - size distribution curve, and ω1, ω2, ω3 are the weights determined by the entropy - weight method.

[0043] Furthermore, for the finite - element numerical model after assignment, adopt an adaptive mesh refinement technique to improve the mesh resolution of the finite - element numerical model and obtain a numerical model with a refined mesh, including:

[0044] In the finite element numerical model after assignment, adaptive mesh refinement is adopted. The global basic mesh size is 1m×1 m×1 m, and mesh subdivision is dynamically triggered according to the following conditions:

[0045] If the displacement gradient threshold ≥5mm / m or the safety factor sensitivity δF s ≥0.01, then mesh subdivision is triggered;

[0046] Each time the mesh is subdivided, the mesh size is halved until the local mesh size is optimized to the 0.05m level, and a numerical model with a refined mesh is obtained.

[0047] Furthermore, the displacement gradient threshold is calculated by the formula:

[0048]

[0049] where , and respectively represent the three components of the displacement field in the x, y, and z directions in three-dimensional space, , and respectively represent their rates of change in the x, y, and z directions.

[0050] The safety factor sensitivity δF s is calculated by the formula:

[0051]

[0052] where pi is the parameter to be analyzed, Fs is the safety factor under the current parameter, represents the rate of change of the safety factor relative to the parameter to be analyzed.

[0053] Furthermore, the method further includes:

[0054] By introducing a GPU-accelerated parallel octree topology algorithm, dynamic partition management of the computational domain is realized. The GPU-accelerated parallel octree topology algorithm and the adaptive mesh refinement method have a cooperative relationship between the underlying data structure and the upper-layer calculation strategy. The nodes (cubes) of the octree correspond to the mesh elements in the computational domain. By dynamically subdividing or coarsening the nodes, adaptive encryption / sparsity of the mesh is directly achieved; the GPU thread pool is used to synchronously process node splitting / merging to accelerate topology update. The octree nodes are divided into blocks according to their spatial positions and assigned to different GPU stream processors (SMs) to maximize parallelism.

[0055] Furthermore, by introducing a GPU-accelerated parallel octree topology algorithm, dynamic partition management of the computational domain is realized, including:

[0056] Build a linear octree data structure, and use Morton coding (Z-order curve) to store the space index. The coding formula is:

[0057]

[0058] where x i , y i , z i are the binary bits of the coordinates of the i-th digit.

[0059] Use CUDA Unified Memory to divide three memory pools: Node pool: Store the octree node attributes, including coordinates, dimensions, parent / child pointers; Field variable pool: Store physical quantities such as displacement d, stress σ, and safety factor Fs; Task pool: Dynamically store the encryption task queue of the mesh to be refined.

[0060] Furthermore, according to the numerical model after refining the mesh, based on the improved NSGA-II algorithm, synchronously generate two groups of candidate slip surfaces, and fuse the strength reduction method and the energy method to construct a collaborative safety factor. Search for parallel double slip surfaces in the two groups of candidate slip surfaces, including:

[0061] According to the refined network, based on the improved NSGA-II algorithm, synchronously generate two groups of candidate slip surfaces, including slip surface A and slip surface B;

[0062] Calculate the safety factor F A of slip surface A and the safety factor F B of slip surface B respectively, and fuse the strength reduction method and the energy method to construct a collaborative safety factor.

[0063] Respectively take the minimum collaborative safety factor and the maximum sliding volume as the objective functions to ensure coverage of the potential dangerous area, and search for parallel double slip surfaces.

[0064] Furthermore, randomly generate 50% arc-shaped and 50% non-arc-shaped slip surface candidate solutions, where:

[0065] For the generation of arc parameters: The abscissa of the center of the circle is randomly selected within the range from the toe of the slope to 2 times the slope height behind the crest of the slope;

[0066] The ordinate of the center of the circle is randomly generated between 30% of the slope height and 1.5 times the slope height; The radius R satisfies the geological rationality constraint;

[0067] For the generation of non-arc parameters: Uniformly distribute 5 control points along the transverse direction of the slope; The ordinates of the control points are randomly perturbed within the range of 20% of the slope height above and below the slope line;

[0068] Convert the parameters into a chromosome vector with real - number encoding to form the initial population. Then, iteratively reduce the soil strength parameters \(c\) and \(\varphi\) according to the geometric parameters of a single sliding surface until the slope displacement mutates. The criterion is: the maximum nodal displacement increment \(> 5\%\) of the slope height, and record the critical reduction factor \(F\). SRM Calculate the cumulative plastic strain energy \(U\) within the sliding zone. plastic Calculate the total strain energy \(U\). total \(=U\). elastic \(+U\). plastic Calculate the energy index \(F\). energy \(=1 - U\). plastic / \(U\). total And calculate the composite coefficient \(F\). composite :

[0069]

[0070] where \(\sigma\). c is the coefficient of variation of the strength parameter, and \(\sigma\). E is the dispersion of the energy distribution.

[0071] Calculate the earthwork volume \(V\) by integrating the area enclosed by the sliding surface and the slope surface, and output the result as a two - dimensional objective space \((1 / F\). composite , \(V)\) standardized from the objective function vector \([1 / F\). composite , \(V]\). Generate \(12\times12\) uniform grid points within the space, associate each candidate solution with the nearest reference point, retain at most 2 optimal solutions for each reference point, preferentially select solutions in sparse regions to avoid local aggregation, and perform adaptive genetic operations:

[0072] Geometric crossover (circular arc type):

[0073] Linear interpolation of the parent - generation center - point coordinates: ;

[0074] The radius takes the weighted average: ;

[0075] Control - point crossover (non - circular arc type):

[0076] Randomly select the crossover position in the parent - generation control - point sequence;

[0077] Exchange the coordinates of the second - half control points and perform cubic - spline smoothing;

[0078] And perform the following mutation operations: Curvature - driven mutation: Calculate the curvature \(k\) of each point on the sliding surface. i , and apply Gaussian perturbation to the control points in the high - curvature region (\(k\). i >\(k\). mean \(+ 2\sigma\). k ). Depth - constraint mutation: If the bottom of the sliding surface exceeds the bedrock depth, lift the deepest control point above the bedrock surface.

[0079] After that, the Pareto solution set is screened: retain the potential dangerous solutions with F composite < 1.05; arrange the candidate solutions in descending order of sliding volume; merge the similar solutions with a spacing less than 0.2 times the slope height; select the top 10 high-volume solutions from the Pareto solution set, and form candidate pairs by combining them in pairs, which are respectively denoted as sliding surface A and sliding surface B.

[0080] Furthermore, the safety factor F A of the sliding surface A and the safety factor F B of the sliding surface B are respectively determined by the following formulas:

[0081]

[0082]

[0083] where c and φ are the strength parameters of the sliding soil mass, L is the segmented length of the sliding surface, N and T are the normal force and the tangential force. Subscripts A and B represent the sliding surface A and the sliding surface B respectively, which are the two searched sliding surfaces.

[0084] Furthermore, the collaborative safety factor F AB in step 5 takes the value formula as:

[0085]

[0086] where the weight coefficients:

[0087] (proportion of shallow sliding volume)

[0088] (proportion of deep sliding volume)

[0089] where , represents the interaction factor, Δz is the vertical spacing between the shallow and deep sliding surfaces, H is the slope height, is the interaction intensity index, V1 is the shallow sliding volume, and V2 is the deep sliding volume.

[0090] Furthermore, according to the obtained parallel double sliding surfaces, the failure probability is solved by Monte Carlo simulation, and the double critical sliding surfaces are judged, including:

[0091] Solve the failure probabilities P A and P B of the double sliding surfaces respectively through Monte Carlo simulation, and calculate the collaborative failure probability P AB :

[0092] If the collaborative safety factor F AB < 1.0, and the collaborative failure probability P ABIf it is > 60%, it is determined that the parallel double sliding surface is a double critical sliding surface, and the loop ends. Otherwise, dynamic grid adaptive adjustment is performed, and the adjustment formula is:

[0093] ;

[0094] where h new is the new grid size, h old is the old grid size, β = 0.1 is the adjustment coefficient, is the modulus of the safety factor gradient, and its formula is as follows,

[0095] ;

[0096] In the formula, , , are the local change rates of F AB in the x, y, and z directions.

[0097] Subsequently, according to the adjusted grid, the sliding surface search range is updated, and the double sliding surface parameter combination is continuously updated until the safety factor F AB converges to a stable minimum value, the sliding surface with the minimum global safety factor is determined, and F A and F B the optimal solutions of the two critical sliding surfaces are locked to obtain a new parallel double sliding surface;

[0098] The failure probabilities P A and P B of the new parallel double sliding surfaces are re-solved. If the cooperative safety factor F AB < 1.0 and the cooperative failure probability P AB > 60%, it is determined that the parallel double sliding surface is a double critical sliding surface. Otherwise, jump back to perform dynamic grid adaptive adjustment again.

[0099] Furthermore, the value formulas of the failure probabilities P A , P B are:

[0100]

[0101]

[0102] The value formula of the said cooperative failure probability P AB is:

[0103]

[0104] where N is the number of sampling times. Here, N represents the total number of Monte Carlo simulations, and I(.) is the indicator function (1 when the condition is satisfied). F A composite,iDenote the composite coefficient of the \(i\)-th Monte Carlo simulation in the sliding surface A, \(F\) B composite,i Denote the composite coefficient of the \(i\)-th Monte Carlo simulation in the sliding surface B.

[0105] Furthermore, the method for judging that the safety factor \(F\) AB converges to a stable minimum value includes:[[]]

[0106] By recording the \(F\) value in each iteration process, when the relative change rate of \(F\) AB ≤ 0.5% and the absolute change ≤ 0.01, it is considered to converge to a stable minimum value. AB

[0107] Furthermore, based on the double critical sliding surface, using the collaborative framework of particle swarm optimization and tabu search to lock the optimal critical sliding surface includes:[[]]

[0108] Using the particle swarm optimization algorithm and the tabu search algorithm to alternately dominate the process of searching for the sliding surface.[[]]

[0109] The alternating trigger conditions include:[[]]

[0110] Periodic switching: Switch the dominant algorithm after every fixed number of iterations.[[]]

[0111] Convergence state trigger: When it is detected that the optimal solution in a continuous preset number of generations has not been improved, automatically switch to another algorithm.[[]]

[0112] In the global exploration stage of the particle swarm: Randomly generate \(N\) particles (such as 50), including circular arc and non-circular arc types, and the velocity update adopts:[[]]

[0113]

[0114] Among them,[[]] represents the updated velocity of particle \(i\) in the \(d\)-th dimension,[[]] is the current position of particle \(i\) in the \(d\)-th dimension,[[]] , , are acceleration coefficients,[[]] =[[]] = 1.5,[[]] = 0.5,[[]] , , are random factors,[[]] , , ~U(0, 1), \(\omega(t)\) is the dynamic inertia weight (initially 0.9, linearly decreasing to 0.4),[[]] is the individual historical best position of particle \(i\) in the \(d\)-th dimension,[[]] is the optimal solution among all particles in the \(d\)-th dimension, \(t\)[[]] best ​It is the current optimal solution passed to the TS stage.

[0115] The speed update adopts:

[0116]

[0117] The fitness evaluation adopts:

[0118]

[0119] Among them, F composite is the collaborative safety factor (strength reduction method + energy method), V is the sliding volume, which needs to be maximized (converted to a minimization problem by taking the reciprocal), and α is the dynamic weight (initially 0.8, decreasing to 0.5 with iteration).

[0120] In the Tabu Search (TS) stage: input the global optimal solution g of PSO best ; generate neighborhood solutions: arc type: perturb the center of the circle (Δx c ∈[−0.1H, 0.1H]) and the radius (ΔR = ±0.1R), and non-arc type: apply Gaussian perturbation (σ = 0.05H) to the point with the maximum curvature; among them, the tabu object is the sliding surface type + key parameter hash (such as rounding the center coordinates to 0.1H); the tabu length is dynamically adjusted (5 - 20 iterations), if the candidate solution is better than the historical optimum, lift the tabu, generate 20 neighborhood solutions in each round, and select the non-tabu solutions to update preferentially.

[0121] Establish an information sharing mechanism between Particle Swarm Optimization (PSO) and Tabu Search (TS). When PSO → TS, inject the top 10% of the elite solutions of the PSO population into the TS candidate pool every 5 generations; when TS → PSO, replace the worst 10% of the particles in PSO with the optimization results of TS, and the tabu list feedback restricts repeated searches.

[0122] The convergence criterion is: the improvement amplitude of the global optimal solution in 20 consecutive generations < 0.1%; the total number of iterations ≥ 200 or the calculation time > 6 hours, lock the optimal double-critical sliding surface morphology and position.

[0123] Beneficial effects: The present invention has the following advantages compared with the prior art:

[0124] (1) By constructing a collaborative search framework for double-critical sliding surfaces, this method innovatively integrates multi-objective optimization algorithms, dynamic grid subdivision technology, and physical mechanism constraints to quantitatively characterize the interaction effect of sliding surfaces and efficiently capture the global optimal solution. This method not only breaks through the mechanical simplification shackles of traditional single-sliding surface models, but also significantly improves the calculation efficiency and interpretation reliability by introducing a data-physics hybrid driving strategy, providing a theoretical breakthrough for the multi-scale analysis of complex slope instability mechanisms and disaster prevention and control, and having important value for promoting the paradigm transformation of geotechnical engineering safety assessment from experience-driven to mechanism-driven.

[0125] (2)Through multi-factor parametric processing and flexible sliding surface morphology adaptation, the present invention can seamlessly connect to various geological terrains and load conditions. From soft clay slopes to high-steep hard rock slopes, and from static loads to seismic dynamic loads, it can accurately capture the double critical state, with a wide range of application scenarios.

[0126] (3)Through fine geological modeling and collaborative safety factor calculation, the present invention accurately depicts the complex mechanical behavior of slopes. Compared with traditional methods, it reduces the position error of the double critical sliding surface and improves the accuracy of the safety factor, providing a reliable quantitative basis for engineering.

[0127] (4)Through intelligent initialization and efficient collaborative search mechanism, the present invention significantly reduces computational redundancy. Facing large-scale slope engineering, the calculation time is shortened and the efficiency is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0128] Figure 1 is the flowchart of the method of the present invention;

[0129] Figure 2 is the schematic diagram of the double critical sliding surface discrimination criterion based on the safety factor and volume;

[0130] Figure 3 is the schematic diagram of the geometric model and soil parameters of the double-layer slope;

[0131] Figure 4 is the double critical sliding surface of the slope in Embodiment 1;

[0132] Figure 5 is the schematic diagram of the three-dimensional slope model;

[0133] Figure 6 is the schematic diagram of the three-dimensional slope mesh division;

[0134] Figure 7 is the schematic diagram of the three-dimensional slope critical sliding surface. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0135] The present invention will be further described below with reference to the drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and cannot be used to limit the protection scope of the present invention.

[0136] As Figure 1 shown is the flowchart of the method of the present invention. It can be seen that a search method for determining the double critical sliding surface of a slope according to the present invention comprises the following steps: Figure 1 It can be seen that a search method for determining the double critical sliding surface of a slope according to the present invention comprises the following steps:

[0137] (1)Obtain geological parameters and topographic data: Through on-site geological exploration, obtain the geotechnical stratification data (cohesion c, internal friction angle φ, unit weight γ), the spatial distribution of weak interlayers (thickness, dip angle), and the dynamic data of the groundwater level (depth, permeability coefficient k) in the study area. Use the vane shear test (for soft soil) and point load test (for rock mass) to correct the errors of laboratory parameters. Use UAV photogrammetry and airborne LiDAR to obtain the digital orthophoto map, digital elevation model (DEM), and digital surface model (DSM) of the study area;

[0138] Among them, the aerial remote sensing digital results obtained by UAV photogrammetry and airborne LiDAR include: Digital Elevation Model (DEM, resolution 0.1m), Digital Surface Model (DSM, resolution 0.1m), Digital Orthophoto Map (DOM, resolution 0.05m), and real-scene three-dimensional model (centimeter level). Through three-dimensional stereoscopic measurement and elevation difference model, accurately extract the geometric shapes and distribution characteristics such as the length, width, area, and volume of the landslide body.

[0139] (2)Construct a three-dimensional real-scene model with the obtained geological parameters and topographic data. Import the constructed model into finite element software to generate a finite element numerical model. Conduct mesh division on the finite element numerical model to transform the continuous solid model into a discrete model composed of a finite number of elements and nodes.

[0140] (3)Set the value range and probability distribution function for the geological parameters obtained in step 1. Use Monte Carlo random simulation to generate an initial parameter sample set with N = 10 4 At the same time, embed the genetic algorithm gene coding, map the parameter group to chromosome individuals, avoid the local convergence caused by the traditional fixed parameter starting point, and randomly select M = 500 groups of parameter combinations from them to ensure that the samples cover the entire domain;

[0141] Among them, the set target parameters are cohesion c (kPa), internal friction angle φ (°), unit weight γ (kN / m³), and pore water pressure u (kPa).

[0142] Table 1 Value range table of geological parameters

[0143]

[0144] For each parameter, a random number is generated according to the set probability distribution function. The random number generation function in the programming language can be used to synthesize the generated random array into a parameter group (in the genetic algorithm, each parameter combination is encoded as a chromosome individual, and the "parameter group" represents a parameter sample) to form an initial parameter sample set. Each parameter group represents a possible parameter combination.

[0145] Determine the encoding method of chromosomes. Map the parameter group to individual chromosomes, encode each parameter group (c, φ, γ, u) as a binary chromosome, and set the gene length according to the parameter accuracy requirements:

[0146] Table 2 Chromosome coding table

[0147]

[0148] Chromosome structure: total length 8+6+4+7=25 bits

[0149] From 10 4 Randomly select M=500 groups from the Monte Carlo samples to ensure that all regions of the parameter space are covered and samples that do not meet the physical constraints (such as c<15kPa or φ>40°) are removed. Then Latin hypercube sampling (LHS) is used to stratify the selected samples to ensure uniform distribution of each parameter dimension. Then, 10% of the individuals are randomly flipped by 1-2 bits to introduce small-scale mutations. The initial population size P=200 is formed, and the chromosome set Population={Chromosome1, Chromosome2, ..., Chromosome 200}.

[0150] 500 groups were randomly selected from the 104 groups for geological matching evaluation. This number (500) is an artificially set fixed value, aiming to balance the contradiction between computational efficiency and parameter space coverage.

[0151] Construct a multi-metric similarity function:

[0152]

[0153] The optimal parameter group is screened out, and different rock and soil layers are assigned values ​​in the numerical model; and the finite element numerical model after the assignment is obtained.

[0154] (4) In the finite element numerical model after assignment, an adaptive grid is used to refine the global basic grid. The initial global grid is divided into 1×1×1m cubic grids, a linear octree data structure is established, and Morton coding (Z-order curve) is used to store the spatial index. The coding formula is:

[0155]

[0156] where x i , y i , z i are the binary bits of coordinates.

[0157] Use CUDA Unified Memory to divide three memory pools: Node pool: Store the attributes of octree nodes (coordinates, sizes, parent / child pointers); Field variable pool: Store physical quantities such as displacement d, stress σ, and safety factor Fs; Task pool: Dynamically store the encryption task queue of the mesh to be refined. Then perform the calculation of the adaptive trigger condition:

[0158] where the displacement gradient threshold ( ) is calculated by the formula:

[0159]

[0160] The sensitivity of the safety factor (δF s ) is calculated by the formula:

[0161]

[0162] where pi is the parameter to be analyzed (such as cohesion c, internal friction angle, etc.), and Fs is the safety factor under the current parameter.

[0163] Step 4 optimizes the local mesh resolution of the finite element model through adaptive mesh refinement (for example, the mesh size near the potential slip surface is refined to the 0.05m level). The generated high-precision mechanical field data (such as displacement gradient, safety factor sensitivity) is more conducive to the search for the candidate set of slip surfaces in Step 5, making the position of the searched double slip surface more accurate.

[0164] This step performs mesh division on the finite element numerical model, and its results are:

[0165] 1. Locally refined mesh: The mesh resolution in the key area reaches the 0.05m level, ensuring the calculation accuracy of mechanical responses;

[0166] Through the adaptive subdivision technology, the basic mesh is optimized from the global 1m×1m×1m to the local 0.05m level, significantly improving the mesh resolution in the key areas (such as areas with large displacement gradients and strong safety factor sensitivities).

[0167] The refined mesh provides high-precision mechanical field data (such as displacement, stress) for Step 5, ensuring that the improved NSGA-II algorithm can accurately calculate the safety factor (F A , F B ) and the collaborative safety factor (F AB ).

[0168] 2. Parallel computing framework: Accelerate octree management through GPU to achieve efficient dynamic partitioning and resource allocation;

[0169] 3. Pre-computed mechanical field data: Displacement gradient and sensitivity field distribution to guide subsequent slip surface search.

[0170] If the displacement gradient threshold ( ) ≥ 5 mm / m and the safety factor sensitivity (δF s ) ≥ 0.01, encryption is triggered, and the local grid size is optimized to the 0.05 m level. Complete the parallel octree topology algorithm accelerated by GPU to achieve dynamic partitioning management of the computational domain;

[0171] (5) Based on the refined grid in step 4, two sets of slip surface candidate sets (slip surface A, slip surface B) are generated synchronously based on the improved NSGA-II algorithm. NSGA-II is the second generation of the Non-dominated Sorting Genetic Algorithm, a classic algorithm for solving multi-objective optimization problems. By randomly generating 50% circular arc type and 50% non-circular arc type slip surface candidate solutions, where:

[0172] For the generation of circular arc parameters: The abscissa of the center of the circle is randomly selected within the range from the toe of the slope to twice the slope height behind the crest;

[0173] The ordinate of the center of the circle is randomly generated between 30% of the slope height and 1.5 times the slope height; the radius R satisfies the geological rationality constraint;

[0174] For the generation of non-circular arc parameters: Five control points are evenly distributed horizontally along the slope surface; the ordinates of the control points are randomly perturbed within the range of 20% of the slope height above and below the slope surface line;

[0175] The parameters are transformed into a chromosome vector with real number encoding to form an initial population. Then, according to the geometric parameters of a single slip surface, the soil strength parameters c and φ are iteratively reduced until the slope displacement mutates (criterion: the maximum node displacement increment > 5% of the slope height), and the critical reduction coefficient F SRM is recorded; calculate the cumulative plastic strain energy U plastic within the slip zone, calculate the total strain energy U total = U elastic + U plastic , calculate the energy index F energy = 1 - U plastic / U total , and perform compounding to obtain the compound index F composite :

[0176]

[0177] where σc is the coefficient of variation of the strength parameter, σ E is the energy distribution dispersion.

[0178] Calculate the earthwork volume V by integrating the area enclosed by the integral sliding surface and the slope surface, and output the result as the objective function vector [1 / F composite , V] Generate 12×12 uniform grid points in the standardized two-dimensional objective space (1 / F composite , V), associate each candidate solution with the nearest reference point, retain at most 2 optimal solutions for each reference point, preferentially select solutions in the sparse area to avoid local aggregation, and perform adaptive genetic (adaptive crossover) operations:

[0179] Geometric crossover (arc type):

[0180] Linear interpolation of the center coordinates of the parent circles: ;

[0181] The radius is taken as the weighted average: ;

[0182] Control point crossover (non-arc type):

[0183] Randomly select the crossover position in the sequence of control points of the parent;

[0184] Exchange the coordinates of the latter half of the control points and perform cubic spline smoothing;

[0185] And perform the following mutation operations: Curvature-driven mutation: (Calculate the curvature k of each point on the sliding surface i , and apply Gaussian perturbation (σ = 0.1H) to the control points in the high-curvature area (k i >k mean +2σ k )); Depth constraint mutation: (If the bottom of the sliding surface exceeds the bedrock depth, lift the deepest control point above the bedrock surface).

[0186] Then perform Pareto solution set screening: Retain potential dangerous solutions with F composite <1.05; Arrange the candidate solutions in descending order of sliding volume; Merge similar solutions with a distance less than 0.2 times the slope height; Select the top N (usually N = 10) high-volume solutions from the Pareto solution set, and form candidate pairs by combining them in pairs (denoted as sliding surface A and sliding surface B respectively), and the corresponding safety factors are F A and F B ;

[0187]

[0188]

[0189] Among them, c and φ are the strength parameters of the sliding soil mass, L is the segmented length of the sliding surface, and N and T are the normal force and the tangential force.

[0190] And calculate the collaborative safety factor F AB :

[0191]

[0192] Among them, the weight coefficient

[0193] (proportion of the shallow sliding volume)

[0194] (proportion of the deep sliding volume)

[0195] (interaction factor, vertical distance Δz between the shallow and deep sliding surfaces, and H is the slope height)

[0196] Respectively, with the minimum safety factor and the maximum sliding volume as the objective functions, to ensure covering the potential dangerous area, search for the parallel double sliding surfaces, as Figure 2 shown. In the figure, calculate the double critical sliding surfaces based on the safety factor and volume discrimination for two typical cases; the potential dangerous area refers to all possible unstable ranges in the landslide body, that is, the potential sliding mass, which is generally judged according to the slope.

[0197] The calculation formula for the sliding volume is:

[0198]

[0199] Among them, A i is the cross-sectional area of the sliding body, and Δz i is the longitudinal grid spacing.

[0200] (6) Solve the failure probabilities P A and P B ,

[0201]

[0202]

[0203] And calculate the collaborative failure probability P AB :

[0204]

[0205] Among them, N is the number of sampling times, I() is the indicator function, which is 1 when the condition is satisfied, and F A composite,i represents the composite coefficient of the i-th Monte Carlo simulation in the sliding surface A, and F Bcomposite,i Represents the composite coefficient of the i-th combination in the sliding surface B.

[0206] If the collaborative safety factor (F AB < 1.0), and the collaborative failure probability P AB > 60%, it is determined as a double-critical sliding surface; otherwise, return to step 4 for dynamic grid adaptive adjustment, and the adjustment formula is:

[0207]

[0208] β = 0.1 is the adjustment coefficient, is the modulus of the safety factor gradient.

[0209] Subsequently, perform step 6 to update the sliding surface search range. Continuously update the parameter combinations of the double sliding surface until the safety factor F AB converges to a stable minimum value, determine the global optimal sliding surface, and lock the optimal solution of the double-critical sliding surface;

[0210] (7) Adopt the collaborative framework of Particle Swarm Optimization (PSO) and Tabu Search (TS):

[0211] In step 7, balance the search efficiency and accuracy by dynamically switching the dominant algorithm and combining the global exploration ability of the Particle Swarm (PSO) and the local optimization ability of the Tabu Search (TS).

[0212] There are two alternating trigger conditions:

[0213] Periodic switching: Switch the dominant algorithm after every fixed number of iterations (such as every 5 generations).

[0214] Convergence state trigger: When it is detected that the current algorithm is stagnant (such as the optimal solution has not improved for 3 consecutive generations), automatically switch to another algorithm.

[0215] In the global exploration stage of the particle swarm: Randomly generate N particles (such as 50), including circular arc and non-circular arc types, and the velocity update adopts:

[0216]

[0217] where ω(t) is the dynamic inertia weight (initially 0.9, linearly decreasing to 0.4), and t best is the current optimal solution passed in the TS stage.

[0218] The velocity update adopts:

[0219]

[0220] The fitness evaluation adopts:

[0221]

[0222] Among them, F composite is the collaborative safety factor (strength reduction method + energy method), V is the sliding volume, which needs to be maximized (converted to a minimization problem by taking the reciprocal), and α is the dynamic weight (initially 0.8, decreasing to 0.5 with iteration).

[0223] In the tabu search (TS) stage: Input the global optimal solution g of PSO best ; Generate neighborhood solutions: Circular arc type: Perturb the center of the circle ( ), and the radius (ΔR = ±0.1R), and non-circular arc type: Apply Gaussian perturbation (σ = 0.05H) to the point with the maximum curvature; among them, the tabu object is the sliding surface type + key parameter hash (such as rounding the center coordinates to 0.1H); the tabu length is dynamically adjusted (5 - 20 iterations), if the candidate solution is better than the historical optimum, lift the tabu, generate 20 neighborhood solutions per round, and select and update the non-tabu solutions preferentially.

[0224] Establish an information sharing mechanism between particle swarm optimization (PSO) and tabu search (TS). When PSO → TS, inject the top 10% of the elite solutions in the PSO population into the TS candidate pool every 5 generations; when TS → PSO, replace the worst 10% of the particles in PSO with the optimization results of TS, and the tabu list feedback restricts repeated searches.

[0225] The convergence criterion is (the improvement amplitude of the global optimal solution in 20 consecutive generations < 0.1%; the total number of iterations ≥ 200 or the calculation time > 6 hours), and lock the optimal double-critical sliding surface form and position.

[0226] This method provides a theoretical basis for the stability evaluation of heterogeneous slopes under complex geological conditions, early detection of potential slope instability risks and landslide disaster warnings, guides engineering design, optimizes design schemes, and reduces engineering costs.

[0227] Example 1

[0228] This example is a heterogeneous double-layer undrained clay slope, and its model is as Figure 3 shown. The undrained shear strength of the upper layer = 60.0 kPa; the internal friction angle = 0°; the unit weight = 20.0 kN / m3; the elastic modulus E = 1×105 Pa; the Poisson's ratio = 0.3; the slope height H = 18.0 m, and these parameters of the slope remain constant. The foundation depth ratio of the slope is still represented by D to indicate the change in foundation depth. The slope model is divided into upper and lower regions. The upper layer is the embankment layer and the lower layer is the foundation layer. The undrained strength of the embankment layer is set to , with a height of H, remains constant. The undrained strength of the foundation layer is assumed to be , and its strength varies and gradually increases. The foundation depth is determined by the variable foundation depth ratio D.

[0229] The non-uniformity of the soil layer strength is represented by the ratio of the foundation layer strength to the undrained strength of the dam layer . The obtained critical slip surface is as shown in Figure 4 , which shows the double critical slip surface of the slope when D = 2.0 and = 2.0, and the strength ratio of the upper and lower soil layers = 1.5.

[0230] The foundation depth ratio D is defined as the relationship ratio between the slope height H and the total slope height (or the foundation layer depth H1):

[0231] Or

[0232]

[0233] Example 2

[0234] This example is a double-layer heterogeneous slope under drained conditions, and its model is as shown in Figure 5 . The strength parameters of the dam layer of the slope are fixed as = 15.0 kPa and = 0.364 ( = 20°). The strength parameters and ( ) of the second foundation layer vary. The ratio of the strength parameters of the two layers is used to measure the strength change of the upper and lower layers. In this example , = 1.19, D = 2.2. In the method described in this article, the three-dimensional slope model is adaptively meshed, as shown in Figure 6 . Then, according to the generated critical slip surface A and slip surface B, their safety factors F A and F B , as well as their cooperative safety factor F AB and cooperative failure probability P AB are calculated. Finally, the position of the double critical slip surface is determined as shown in Figure 7 .

[0235] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0236] The present application is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to the embodiments of the present application. It should be understood that each flow and / or block in the flowchart and / or block diagram can be implemented by computer program instructions, and the combination of flows and / or blocks in the flowchart and / or block diagram can also be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in one Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0237] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including instruction means, and the instruction means implements the functions specified in one Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0238] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process. Therefore, the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0239] The above is only the preferred embodiment of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the technical principle of the present invention, several improvements and modifications can be made, and these improvements and modifications should also be regarded as the protection scope of the present invention.

Claims

1. A search method for determining double critical sliding surfaces of a slope in a heterogeneous soil layer, characterized in that: The following steps are involved: Step 1: Obtain geological parameters and topographic data of the area to be analyzed; Step 2: construct a three-dimensional real scene model and generate a finite element numerical model based on the geological parameters and terrain data; Step 3: Set the parameter value range for the acquired geological parameters, determine the probability distribution function of the parameters, generate and screen parameter samples using Monte Carlo random simulation and genetic algorithm, and determine the optimal parameter group through multi-index similarity function and assign it to the finite element numerical model; Step 4: Adopting adaptive mesh refinement technology to improve the mesh resolution of the finite element numerical model after the assignment, and obtaining the numerical model after mesh refinement; Step 5: According to the numerical model after mesh refinement, two sets of sliding surface candidate sets are generated synchronously based on the improved NSGA-II algorithm, and the strength reduction method and energy method are integrated to construct the collaborative safety factor. The parallel dual sliding surfaces are searched in the two sets of sliding surface candidate sets; Step 6: Based on the obtained parallel dual sliding surfaces, the failure probability is solved by Monte Carlo simulation, and the dual critical sliding surfaces are determined; Step 7: Based on the dual critical sliding surface, the particle swarm optimization and taboo search collaborative framework is used to lock the optimal critical sliding surface.

2. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: Obtain geological parameters and topographic data, including: Investigate geological parameters of the geology through on-site geological exploration, including rock and soil stratification parameters, spatial distribution of weak interlayers and groundwater level dynamic data; The rock and soil stratification parameters include cohesion c, internal friction angle φ, and bulk density γ; The spatial distribution of the weak interlayer includes the thickness and inclination of the weak interlayer space; The groundwater level dynamic data include groundwater depth and permeability coefficient k; UAV photogrammetry and airborne LiDAR were used to obtain terrain data, including digital orthophotos, digital terrain models (DEMs), and digital surface models (DSMs) of the study area.

3. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: According to the geological parameters and terrain data, a three-dimensional real scene model is constructed and a finite element numerical model is generated, including: Construct a three-dimensional real scene model using the acquired terrain data; Import the constructed 3D real-scene model into the finite element software to generate a finite element numerical model. The initial meshing is performed on the entire area of ​​the finite element numerical model, and the entire finite element numerical model is discretized into a finite number of cells and nodes with equal meshes.

4. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: The parameter value range is set for the acquired geological parameters, the probability distribution function of the parameters is determined, the parameter samples are generated and screened using Monte Carlo random simulation and genetic algorithm, and the optimal parameter group is determined through the multi-index similarity function and assigned to the finite element numerical model, including: The parameter value range and probability distribution function are set for the acquired geological parameters, and the parameter groups are generated by Monte Carlo random simulation to form the initial parameter sample set. Each parameter group represents a parameter combination, and the genetic algorithm gene coding is embedded to map the parameter group to chromosome individuals. Randomly select M groups of parameters from the initial parameter sample set to ensure that the samples cover the entire domain of the finite element numerical model. The geological matching degree of each selected parameter combination is evaluated through a multi-index similarity function to select the optimal parameter group; Different geotechnical layers are assigned values ​​in the finite element numerical model according to the optimal parameter set.

5. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 4, characterized in that: The multi-index similarity function is: ； in, , represents the normalized cohesion difference, where represents the sample cohesion, Represents the actual measured cohesion on site, represents the standard deviation of the field cohesion, is the internal friction angle deviation angle, KS(GSD) is the Kolmogorov-Smirnov test statistic of the particle grading curve, ω1, ω2, ω3 are the weights determined by the entropy weight method.

6. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: For the finite element numerical model after assignment, the adaptive mesh refinement technology is used to improve the mesh resolution of the finite element numerical model, and the numerical model after mesh refinement is obtained, including: Adaptive mesh refinement is used in the assigned finite element numerical model, and mesh subdivision is dynamically triggered according to the following conditions: If the displacement gradient threshold ≥5mm / m or safety factor sensitivity δF s ≥0.01 triggers mesh subdivision; Each time the mesh is subdivided, the mesh size is halved until the local mesh size is optimized to the 0.05m level, and the numerical model after the mesh is refined is obtained; Displacement gradient threshold The value formula is: ； in , and They represent the three components of the displacement field along the x, y, and z directions in three-dimensional space, , and Respectively represent the rate of change along the x, y, and z directions; Safety factor sensitivity δF s The value formula is: ； Where pi is the parameter to be analyzed, Fs is the safety factor under the current parameters, It represents the rate of change of the safety factor relative to the parameter to be analyzed.

7. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: According to the numerical model after mesh refinement, two sets of sliding surface candidate sets are generated synchronously based on the improved NSGA-II algorithm, and the synergistic safety factor is constructed by integrating the strength reduction method and the energy method. The parallel dual sliding surfaces are searched in the two sets of sliding surface candidate sets, including: According to the refined network, based on the improved NSGA-II algorithm, two sets of sliding surface candidate sets are synchronously generated, including sliding surface A and sliding surface B; Calculate the safety factor F of sliding surface A separately A and safety factor F of sliding surface B B , integrating the strength reduction method and the energy method to construct a synergistic safety factor, The minimum cooperative safety factor and the maximum sliding volume are taken as objective functions to ensure that the potential dangerous area is covered, and the parallel double sliding surfaces are searched.

8. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 7, characterized in that: The safety factor F of the sliding surface A A and safety factor F of sliding surface B B Determined by the following formulas: ； ； Among them, c, φ are the strength parameters of sliding soil, L is the segment length of sliding surface, N, T are the normal force and tangential force; subscripts A and B represent sliding surface A and sliding surface B respectively, which are the two searched sliding surfaces; Cooperative safety factor F AB The value formula is: ； Among them, the weight coefficient is: ； ； in, , represents the interaction factor, Δz is the vertical distance between the shallow and deep sliding surfaces, H is the slope height, is the interaction strength index, V1 is the shallow sliding volume, and V2 is the deep sliding volume.

9. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: Based on the obtained parallel dual sliding surfaces, the failure probability is solved by Monte Carlo simulation, and the dual critical sliding surfaces are determined, including: The failure probability P of the sliding surface A of each of the two sliding surfaces is solved by Monte Carlo simulation A and the failure probability P of sliding surface B B , and calculate the probability of coordinated failure P AB : If the coordinated safety factor F AB <1.0, and the probability of coordinated failure P AB >60%, the parallel dual sliding surface is determined to be a dual critical sliding surface and the cycle ends. Otherwise, dynamic mesh adaptive adjustment is performed. The adjustment formula is: ； Among them, h new is the new grid size, h old is the old grid size, β is the adjustment coefficient, is the safety factor gradient modulus, and its formula is as follows: ； In the formula, , , F AB The local rate of change in the x, y, and z directions; Then, according to the adjusted grid, the sliding surface search range is updated, and the dual sliding surface parameter combination is continuously updated until the safety factor F AB Converge to a stable minimum, determine the sliding surface with the minimum global safety factor, and lock F A and F B The optimal solution of two critical sliding surfaces gives a new parallel dual sliding surface; Resolve the failure probability P of the new parallel double sliding surfaces A and P B , if the coordinated safety factor F AB <1.0, and the probability of coordinated failure P AB >60%, the parallel dual sliding surface is determined to be a dual critical sliding surface, otherwise the process jumps to re-adjust the dynamic mesh adaptively; Failure probability P of sliding surface A A and the failure probability P of sliding surface B B The value formula is: ； ； The coordinated failure probability P AB The value formula is: ； Where N is the number of sampling times, I() is the indicator function, which is 1 when the condition is met, and F A composite,i represents the composite coefficient of the i-th Monte Carlo simulation in the sliding surface A, F B composite,i represents the composite coefficient of the i-th combination in the sliding surface B.

10. The method for searching for double critical sliding surfaces of a slope in a heterogeneous soil layer according to claim 1, characterized in that: Based on the dual critical sliding surface, the particle swarm optimization and taboo search collaborative framework are used to lock the optimal critical sliding surface, including: The particle swarm optimization algorithm and the taboo search algorithm are used to alternately lead the search for the sliding surface process. Alternate trigger conditions include: Periodic switching: switch the dominant algorithm after completing a fixed number of iterations; or, Convergence state trigger: When it is detected that the optimal solution has not been improved for a preset number of consecutive generations, it automatically switches to another algorithm.

Images