Searching method for determining heterogeneous soil layer slope double-critical sliding surface

Abstract

The invention relates to a searching method for determining a heterogeneous soil layer slope double-critical sliding surface. According to the method, Monte Carlo stochastic simulation is combined with a genetic algorithm to generate and screen parameter samples, and an optimal parameter group is determined through a multi-index similarity function and assigned to a model. And two groups of sliding surface candidate sets are synchronously generated based on an improved NSGA-II algorithm, a collaborative safety coefficient is constructed by fusing a strength reduction method and an energy method, and the side slope mechanical response is comprehensively considered. The failure probability is solved through Monte Carlo simulation, and the dual-critical sliding surface is judged. And finally, locking an optimal critical sliding surface by using a particle swarm optimization and tabu search collaborative framework. The algorithm can accurately capture the double-critical sliding surface, is widely suitable for various geological and load working conditions, improves the calculation precision and efficiency, and provides key data support for slope engineering design, safety evaluation, disaster early warning and the like.

Description

Technical Field

[0001] The invention relates to the technical field of slope stability analysis, and in particular to a search method for determining double critical sliding surfaces of a slope in a heterogeneous soil layer. 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 framework of rigid limit equilibrium theory, and the safety factor of potential sliding surfaces is solved by constructing static equilibrium equations. This type of method simplifies the sliding body into rigid body motion and discretizes the sliding surface using the strip method. Its core assumption is the uniformity of the sliding surface morphology and the idealization of the sliding mechanism. However, actual engineering slopes are often controlled by the coupling of multiple factors: the stratum structure shows significant heterogeneity, and the weak interlayer and crack network form a complex stress transfer path; the spatiotemporal variability of groundwater seepage field and rock and soil mechanical parameters leads to a strength weakening effect; the dynamic application and unloading process of external loads induces progressive failure. Under such complex conditions, the slope instability mode often manifests as a composite failure mechanism with the coordinated action of multiple sliding surfaces, and the single sliding surface assumption is difficult to accurately characterize its mechanical essence.

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

[0004] First, the modeling capability of the interaction effect of sliding surfaces is insufficient, and it is impossible to quantify the energy transfer and kinematic correlation between the two sliding surfaces;

[0005] Second, parameter sensitivity analysis is mostly limited to single-factor disturbances, lacking a global assessment of system responses under multi-parameter coordinated variation;

[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 judgment of the critical state of instability.

[0007] With the development of computational mechanics and optimization algorithms, slope stability analysis methods based on numerical simulation have gradually become a research hotspot. Although numerical methods such as finite element method and discrete element method can more accurately describe the stress and strain field of the slope, they still face severe challenges in the search for critical sliding surfaces.

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

[0009] Deterministic algorithms (such as gradient descent method) rely on initial parameter settings and are prone to fall into local optimal solutions and miss the true critical sliding surface;

[0010] Although random algorithms (such as genetic algorithms and simulated annealing) have the potential for global exploration, the computational cost increases exponentially when the dimension of the parameter space explodes. In addition, multi-sliding surface search requires the construction of a high-dimensional optimization model. The traditional single-objective optimization framework is difficult 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 mesh discretization level, the static uniform partitioning strategy cannot adaptively capture the strain concentration area near the sliding surface, resulting in a waste of computing resources and loss of accuracy.

[0011] Although existing studies have attempted to introduce dynamic grid technology, the subdivision triggering mechanism is mostly based on empirical thresholds and lacks a self-consistent association with the mechanical response characteristics, making it difficult to achieve coordinated optimization of computational efficiency and accuracy. The above bottlenecks restrict the ability to interpret complex slope instability mechanisms, and it is urgent to build a new analysis framework that integrates multi-source data drive, intelligent optimization strategy 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 instability mechanism of multiple sliding surface coupling:

[0013] First, 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 dual sliding surfaces on the stability of the system in actual engineering. Such simplified assumptions lead to the model's inability to capture the stress redistribution, energy transfer and kinematic association between sliding surfaces, resulting in systematic deviations in safety factor assessment, especially in heterogeneous slopes with weak interlayers and interlaced faults, where the risk of misjudgment increases significantly.

[0014] Second, although the existing numerical optimization algorithm has shown certain effectiveness in the search for a single sliding surface, its single-objective optimization mechanism is difficult to coordinate the multi-constraint contradictions between the sliding surface geometry, the minimization of the safety factor and the maximization of the sliding volume, resulting in incomplete coverage of the solution space and missing the real critical sliding surface combination. 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 needs of high-precision engineering decision-making. Summary of the invention

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

[0016] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0017] In a first aspect, the present invention provides a search method for determining a double critical sliding surface of a slope in a heterogeneous soil layer, which specifically comprises the following steps:

[0018] Step 1: Obtain 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 based on the geological parameters and terrain data;

[0020] 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;

[0021] 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;

[0022] 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;

[0023] 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;

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

[0025] Furthermore, geological parameters and topographic data are obtained, including:

[0026] 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;

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

[0028] The spatial distribution of the weak interlayer includes thickness and inclination;

[0029] The groundwater level dynamic data include burial depth and permeability coefficient k.

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

[0031] Furthermore, 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:

[0032] Construct a three-dimensional real scene model using the acquired terrain data;

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

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

[0035] Furthermore, 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 a multi-index similarity function and assigned to the finite element numerical model, including:

[0036] 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 possible parameter combination, and the genetic algorithm gene coding is embedded to map the parameter group to chromosome individuals.

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

[0038] The geological matching degree of each selected parameter combination is evaluated through a multi-index similarity function to select the optimal parameter group;

[0039] Different geotechnical layers are assigned values ​​in the finite element numerical model according to the optimal parameter set.

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

[0041]

[0042] 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 The weights determined by the entropy weight method.

[0043] Furthermore, for the finite element numerical model after the assignment, the adaptive mesh refinement technology is used to improve the mesh resolution of the finite element numerical model, and the numerical model after the mesh refinement is obtained, including:

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

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

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

[0047] Furthermore, the displacement gradient threshold The formula for obtaining the value is:

[0048]

[0049] in , and They represent the three components of the displacement field along the x, y, and z directions in three-dimensional space, , and Represent the rate of change along the x, y, and z directions respectively.

[0050] Safety factor sensitivity δF s The formula for obtaining the value is:

[0051]

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

[0053] Furthermore, the method further comprises:

[0054] By introducing the 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 subdivision method are the collaborative relationship between the underlying data structure and the upper-level computing strategy. The nodes (cubes) of the octree correspond to the grid units in the computational domain. By dynamically subdividing or coarsening the nodes, the adaptive encryption / sparseness of the grid is directly realized; the GPU thread pool is used to synchronously process node splitting / merging, accelerate topology updates, and divide the octree nodes into blocks according to spatial positions and assign them to different GPU stream processors (SMs) to maximize parallelism.

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

[0056] A linear octree data structure is established, and Morton encoding (Z-order curve) is used to store spatial indexes. The encoding formula is:

[0057]

[0058] where x i ,y i , z i is the coordinate binary bit of the i-th position.

[0059] Use CUDA Unified Memory to divide the memory into three pools: node pool: stores octree node attributes, including coordinates, size, parent / child pointers; field variable pool: stores physical quantities such as displacement d, stress σ, safety factor Fs, etc.; task pool: dynamically stores the encrypted task queue of the grid to be subdivided.

[0060] Furthermore, 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 the energy method are integrated to construct the synergistic safety factor. The parallel dual sliding surfaces are searched in the two sets of sliding surface candidate sets, including:

[0061] 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;

[0062] 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,

[0063] The minimum cooperative safety factor and the maximum sliding volume are taken as objective functions respectively to ensure the coverage of potential dangerous areas and to obtain the parallel double sliding surfaces through search.

[0064] Furthermore, 50% circular arc and 50% non-circular arc sliding surface candidate solutions are randomly generated, where:

[0065] For arc parameter generation: the horizontal coordinate of the center of the circle is randomly selected within the range of twice the slope height from the slope foot to the slope top;

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

[0067] For non-circular arc parameter generation: 5 control points are evenly distributed along the slope surface; the vertical coordinates of the control points are randomly disturbed within the range of 20% of the slope height above and below the slope line;

[0068] The parameters are converted into chromosome vectors encoded by real numbers to form the initial population. Then, the soil strength parameters c and φ are iteratively reduced according to the geometric parameters of a single sliding surface until the slope displacement suddenly changes. The criterion is: the maximum node displacement increment is greater than 5% of the slope height, and the critical reduction factor F is recorded. SRM ; Calculate the accumulated plastic strain energy U in 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 σ c is the coefficient of variation of the strength parameter, σ E is the energy distribution dispersion.

[0071] The earthwork volume V is calculated by integrating the area enclosed by the sliding surface and the slope surface, and the result is output as the objective function vector [1 / F composite , V] The standardized two-dimensional target space (1 / F composite , V), generate 12×12 uniform grid points, associate each candidate solution to the nearest reference point, retain at most 2 optimal solutions for each reference point, give priority to solutions in sparse areas, avoid local aggregation, and perform adaptive genetic operations:

[0072] Geometric intersection (arc type):

[0073] Linear interpolation of parent circle center coordinates: ;

[0074] The radius is weighted average: ;

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

[0076] The crossover position is randomly selected in the parent control point sequence;

[0077] Exchange the coordinates of the control points in the second half 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 , for high curvature regions (k i >k mean +2σ k ) is subjected to Gaussian perturbation; Depth-constrained variation: if the bottom of the sliding surface exceeds the bedrock burial depth, the deepest control point is lifted above the bedrock surface.

[0079] Then perform Pareto solution screening: retain F composite <1.05 potential dangerous solutions; arrange candidate solutions from large to small according to sliding volume; merge 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, combine them in pairs to form candidate pairs, and record them as sliding surface A and sliding surface B respectively.

[0080] Furthermore, the safety factor F of the sliding surface A is A and safety factor F of sliding surface B B Determined by the following formulas:

[0081]

[0082]

[0083] Among them, c, φ are the strength parameters of sliding soil, L is the length of the sliding surface segment, N, T are the normal force and the tangential force. The subscripts A and B represent sliding surface A and sliding surface B, respectively, which are the two searched sliding surfaces.

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

[0085]

[0086] Among them, the weight coefficient is:

[0087] (Volume ratio of shallow sliding)

[0088] (Deep sliding volume ratio)

[0089] 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, V 1 is the shallow sliding volume, V 2 is the deep sliding volume.

[0090] Furthermore, 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:

[0091] The failure probability P of each of the two sliding surfaces is solved by Monte Carlo simulation A and P B , and calculate the probability of coordinated failure P AB :

[0092] 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:

[0093] ;

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

[0095] ;

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

[0097] 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;

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

[0099] Furthermore, the failure probability P A , P B The value formula is:

[0100]

[0101]

[0102] The coordinated failure probability P AB The value formula is:

[0103]

[0104] Where N is the number of samples, where N represents the total number of Monte Carlo simulations, and I (.) is the indicator function (1 when the condition is met).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 Monte Carlo simulation in the sliding surface B.

[0105] Further, determine the safety factor F AB Methods that converge to a stable minimum include:

[0106] By recording the F AB Value, when F AB When the relative change rate is ≤0.5% and the absolute change is ≤0.01, it is considered to have converged to a stable minimum.

[0107] Furthermore, 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:

[0108] The particle swarm optimization algorithm and the taboo search algorithm are used to alternately lead the search for the sliding surface process.

[0109] Alternate trigger conditions include:

[0110] Periodic switching: The dominant algorithm is switched after a fixed number of iterations.

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

[0112] In the global exploration phase of the particle swarm: N particles (e.g. 50) are randomly generated, including arc and non-arc types, and the speed is updated using:

[0113]

[0114] in, represents the updated velocity of particle i in the dth dimension, The current position of particle i in the dth dimension, , , is the acceleration factor, = =1.5, =0.5, , , is a random factor, , , ~U(0, 1), ω(t) is the dynamic inertia weight (initial 0.9, linearly decreasing to 0.4), is the best individual historical position of particle i in the dth dimension, is the optimal solution in the dth dimension among all particles, t best It is the current optimal solution passed in the TS stage.

[0115] Speed ​​update using:

[0116]

[0117] Fitness evaluation uses:

[0118]

[0119] where F composite is the collaborative safety factor (strength reduction method + energy method), V is the sliding volume, which needs to be maximized (taking the inverse to convert it into a minimization problem), and α is the dynamic weight (initial 0.8, decreasing to 0.5 with iterations).

[0120] In the taboo search (TS) phase: input the global optimal solution g of PSO best ; Generate neighborhood solution: Arc type: Perturb the center of the circle (Δx c ∈[−0.1H, 0.1H]) and radius (ΔR=±0.1R) and non-arc type: Gaussian perturbation (σ=0.05H) is applied to the point with maximum curvature; the taboo object is the sliding surface type + key parameter hash (such as the coordinates of the center of the circle are rounded to 0.1H); the taboo length is dynamically adjusted (5~20 iterations), if the candidate solution is better than the historical optimal, the taboo is lifted, 20 neighborhood solutions are generated in each round, and non-tabu solutions are screened and updated selectively.

[0121] An information sharing mechanism between particle swarm optimization (PSO) and taboo search (TS) is established. When PSO→TS, the top 10% elite solutions of the PSO population are injected into the TS candidate pool every 5 generations. When TS→PSO, the worst 10% particles in PSO are replaced by TS optimization results, and the taboo table feedback is used to limit repeated searches.

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

[0123] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0124] (1) This method constructs a dual-critical sliding surface collaborative search framework, innovatively integrates multi-objective optimization algorithms, dynamic mesh subdivision technology and physical mechanism constraints, and realizes quantitative characterization of sliding surface interaction effects and efficient capture of global optimal solutions. This method not only breaks through the mechanical simplification constraints of the traditional single sliding surface model, but also significantly improves the computational efficiency and interpretation reliability by introducing a data-physics hybrid drive strategy, providing a theoretical breakthrough for the multi-scale analysis and disaster prevention of complex slope instability mechanisms, and is of great value in promoting the paradigm shift of geotechnical engineering safety assessment from experience-driven to mechanism-driven.

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

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

[0127] (4) The present invention significantly reduces computational redundancy through intelligent initialization and efficient collaborative search mechanism. When faced with large-scale slope projects, the computational time is shortened and the efficiency is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0128] Figure 1 is a flow chart of the method of the present invention;

[0129] Figure 2 Schematic diagram of double critical sliding surface judgment criteria based on safety factor and volume;

[0130] Figure 3 The geometric model and soil parameter diagram of the double-layer slope;

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

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

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

[0134] Figure 7 Schematic diagram of the critical sliding surface of a three-dimensional slope. DETAILED DESCRIPTION

[0135] The present invention will be further described below in conjunction with the accompanying 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] like Figure 1 Shown is a flow chart of the method of the present invention, comprising Figure 1 It can be seen that the present invention provides a search method for determining the double critical sliding surface of a slope, and the steps are as follows:

[0137] (1) Obtaining geological parameters and topographic data: Through on-site geological exploration, the rock and soil stratification data (cohesion c, internal friction angle φ, bulk density γ), the spatial distribution of weak interlayers (thickness, inclination) and the dynamic data of groundwater level (burial depth, permeability coefficient k) of the study area are obtained. The cross-plate shear test (soft soil) and point load test (rock mass) are used to correct the laboratory parameter errors. UAV photogrammetry and airborne LiDAR are used to obtain digital orthophotos, digital elevation models (DEMs) and digital surface models (DSMs) of the study area.

[0138] The digital results of aerial remote sensing 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-life three-dimensional model (centimeter level). The geometric shape and distribution characteristics of the landslide, such as length, width, area, volume, etc., are accurately extracted through three-dimensional stereo measurement and elevation difference model.

[0139] (2) The acquired geological parameters and terrain data are used to construct a three-dimensional real-scene model, which is then imported into finite element software to generate a finite element numerical model. The finite element numerical model is meshed to convert the continuous solid model into a discrete model consisting of a finite number of units and nodes.

[0140] (3) Set the parameter value range and probability distribution function for the geological parameters obtained in step 1, and use Monte Carlo random simulation to generate N=10 4 The initial parameter sample set is embedded with the genetic algorithm gene code, and the parameter group is mapped to the chromosome individual to avoid the local convergence caused by the traditional fixed parameter starting point. M=500 groups of parameter combinations are randomly selected from them to ensure that the samples cover the entire domain.

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

[0142] Table 1 Geological parameter value ranges

[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={Chromosome 1 , Chromosome 2 , ..., 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 is the coordinate binary bit.

[0157] Use CUDA Unified Memory to divide the memory into three pools: node pool: store octree node attributes (coordinates, size, parent / child pointers); field variable pool: store physical quantities such as displacement d, stress σ, safety factor Fs; task pool: dynamically store the encrypted task queue of the grid to be subdivided. Then perform adaptive trigger condition calculation:

[0158] The displacement gradient threshold ( ) The value formula is:

[0159]

[0160] Safety factor sensitivity (δF s ) The value formula is:

[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 parameters.

[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 sliding surface is refined to the level of 0.05 m). The generated high-precision mechanical field data (such as displacement gradient and safety factor sensitivity) is more conducive to the search for the sliding surface candidate set in step 5, making the searched dual sliding surface positions more accurate.

[0164] This step is to mesh the finite element numerical model, and its results are:

[0165] 1. Locally refined grid: The grid resolution in key areas reaches 0.05m level to ensure the accuracy of mechanical response calculation;

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

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

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

[0169] 3. Pre-calculate mechanical field data: displacement gradient and sensitivity field distribution to guide subsequent sliding surface search.

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

[0171] (5) Based on the refined grid in step 4, two sets of sliding surface candidate sets (sliding surface A and sliding surface B) are generated synchronously based on the improved NSGA-II algorithm. NSGA-II is the non-dominated sorting genetic algorithm II, which is a classic algorithm for solving multi-objective optimization problems. By randomly generating 50% circular and 50% non-circular sliding surface candidate solutions, the following equations are used:

[0172] For arc parameter generation: the horizontal coordinate of the center of the circle is randomly selected within the range of twice the slope height from the slope foot to the slope top;

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

[0174] For non-circular arc parameter generation: 5 control points are evenly distributed along the slope surface; the vertical coordinates of the control points are randomly disturbed within the range of 20% of the slope height above and below the slope line;

[0175] The parameters are converted into real-number coded chromosome vectors to form the initial population. Then, the soil strength parameters c and φ are iteratively reduced according to the geometric parameters of a single sliding surface until the slope displacement suddenly changes (criterion: the maximum node displacement increment > 5% of the slope height), and the critical reduction factor F is recorded. SRM ; Calculate the accumulated plastic strain energy U in 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 compound them to get the composite index F composite :

[0176]

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

[0178] The earthwork volume V is calculated by integrating the area enclosed by the sliding surface and the slope surface, and the result is output as the objective function vector [1 / F composite , V] The standardized two-dimensional target space (1 / F composite , V) generates a 12×12 uniform grid of points, associates each candidate solution to the nearest reference point, retains at most 2 optimal solutions for each reference point, gives priority to solutions in sparse areas, avoids local aggregation, and performs adaptive genetic (adaptive crossover) operations:

[0179] Geometric intersection (arc type):

[0180] Linear interpolation of parent circle center coordinates: ;

[0181] The radius is weighted average: ;

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

[0183] The crossover position is randomly selected in the parent control point sequence;

[0184] Exchange the coordinates of the control points in the second half 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 , for high curvature regions (k i >k mean +2σ k ) is subjected to Gaussian perturbation (σ=0.1H)); depth constraint variation: (if the bottom of the sliding surface exceeds the bedrock burial depth, the deepest control point is raised above the bedrock surface).

[0186] Then perform Pareto solution screening: retain F composite <1.05 potential dangerous solutions; arrange candidate solutions from large to small according to sliding volume; merge similar solutions with a spacing less than 0.2 times the slope height; select the first N (usually N=10) high volume solutions from the Pareto solution set, and combine them in pairs to form candidate pairs (respectively recorded as sliding surface A and sliding surface B), and the corresponding safety factors are F A and F B ;

[0187]

[0188]

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

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

[0191]

[0192] Among them, the weight coefficient

[0193] (Volume ratio of shallow sliding)

[0194] (Deep sliding volume ratio)

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

[0196] Taking the minimum safety factor and the maximum sliding volume as the objective functions, the potential danger area is covered and the parallel double sliding surfaces are searched, such as Figure 2 As shown in the figure, the double critical sliding surface based on safety factor and volume is calculated for two typical cases; the potential danger area refers to all possible instability ranges in the landslide body, that is, the potential sliding body, which is generally judged based on the slope.

[0197] The sliding volume calculation formula is:

[0198]

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

[0200] (6) Determine the failure probability P of each sliding surface through Monte Carlo simulation A and P B ,

[0201]

[0202]

[0203] And calculate the probability of coordinated failure P AB :

[0204]

[0205] 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, FB composite,i represents the composite coefficient of the i-th combination in the sliding surface B.

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

[0207]

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

[0209] Then proceed to step 6 to update the sliding surface search range. Keep updating the dual sliding surface parameter combination until the safety factor F AB Converge to a stable minimum, determine the global optimal sliding surface, and lock the optimal solution of the double critical sliding surface;

[0210] (7) Using the collaborative framework of particle swarm optimization (PSO) and taboo search (TS):

[0211] Step 7 balances search efficiency and accuracy by dynamically switching the dominant algorithm, combining the global exploration capability of particle swarm (PSO) and the local optimization capability of taboo search (TS).

[0212] There are two types of alternating trigger conditions:

[0213] Periodic switching: The dominant algorithm is switched after a fixed number of iterations (such as every 5 generations).

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

[0215] In the global exploration phase of the particle swarm: N particles (e.g. 50) are randomly generated, including arc and non-arc types, and the speed is updated using:

[0216]

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

[0218] Speed ​​update using:

[0219]

[0220] Fitness evaluation uses:

[0221]

[0222] where F composite is the collaborative safety factor (strength reduction method + energy method), V is the sliding volume, which needs to be maximized (taking the inverse to convert it into a minimization problem), and α is the dynamic weight (initial 0.8, decreasing to 0.5 with iterations).

[0223] In the taboo search (TS) phase: input the global optimal solution g of PSO best ; Generate neighborhood solution: Arc type: Perturb the center of the circle ( ) and radius (ΔR=±0.1R) and non-arc type: apply Gaussian perturbation (σ=0.05H) to the point with maximum curvature; the taboo object is the sliding surface type + key parameter hash (such as the coordinates of the center of the circle are rounded to 0.1H); the taboo length is dynamically adjusted (5~20 iterations), if the candidate solution is better than the historical optimal, the taboo is lifted, 20 neighborhood solutions are generated in each round, and non-tabu solutions are screened and updated selectively.

[0224] An information sharing mechanism between particle swarm optimization (PSO) and taboo search (TS) is established. When PSO→TS, the top 10% elite solutions of the PSO population are injected into the TS candidate pool every 5 generations. When TS→PSO, the worst 10% particles in PSO are replaced by TS optimization results, and the taboo table feedback is used to limit repeated searches.

[0225] The convergence criteria are (the improvement of the global optimal solution for 20 consecutive generations is <0.1%; the total number of iterations is ≥200 or the calculation time is >6 hours), and the optimal double critical sliding surface morphology and position are locked.

[0226] This method provides a theoretical basis for the stability assessment of heterogeneous slopes under complex geological conditions, early detection of potential slope instability risks and landslide disaster warning, guiding engineering design, optimizing design solutions and reducing engineering costs.

[0227] Example 1

[0228] This example is a non-homogeneous double-layer undrained clay slope, and its model is as follows: Figure 3 The undrained shear strength of the upper layer =60.0kPa; internal friction angle =0°; bulk density =20.0kN / m3; elastic modulus E=1×105Pa; Poisson's ratio =0.3; slope height H=18.0m, these parameters of the slope remain constant. The foundation depth ratio of the slope is still represented by D to indicate the change of foundation depth. The slope model is divided into two areas, 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 height H, which remains constant. The undrained strength of the foundation layer is assumed to be , its strength is variable and gradually increases, and the foundation depth is determined by the variable foundation depth ratio D.

[0229] The inhomogeneity of soil strength is measured using the ratio of the foundation strength to the undrained strength of the embankment layer. The obtained critical sliding surface is shown as Figure 4 As shown, it shows that D = 2.0 and =2.0, upper and lower soil layer strength ratio =1.5 is the double critical sliding surface of the slope.

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

[0231] or

[0232]

[0233] Example 2

[0234] This embodiment is a double-layer heterogeneous slope under drainage conditions, and its model is as follows Figure 5 As shown, the strength parameter of the embankment layer of the slope is fixed as =15.0kPa and =0.364( =20°), the strength parameters of the second foundation layer and ( ) is varied, applying two layers of intensity parameter ratio To measure the change in the strength of the upper and lower layers, in this embodiment , =1.19, D=2.2, in the method described in this paper, the three-dimensional slope model is adaptively meshed, such as Figure 6 As shown, the safety factor F of the generated critical sliding surface A and sliding surface B is calculated. A and F B , and their coordinated safety factor F AB and the probability of coordinated failure P AB Finally, the position of the double critical sliding surface is determined as Figure 7 shown.

[0235] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application may adopt 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.) containing computer-usable program codes.

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

[0237] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0238] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0239] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection 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, use Monte Carlo random simulation and genetic algorithm to generate and screen parameter samples, 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 respectively to ensure the coverage of potential dangerous areas and to obtain the parallel double sliding surfaces through search.

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