Loess red-bed landslide risk dynamic deduction method based on digital twinning

Abstract

The present application relates to the technical field of computer-aided simulation, in particular to a loess-red bed landslide risk dynamic deduction method based on digital twinning; the method comprises the following steps: obtaining target landslide terrain, loess-red bed contact surface, rear edge crack, potential sliding boundary and disaster-bearing body data, and constructing a digital twinning contact surface model; calculating the water stagnation height according to the rainfall intensity, the rear edge crack connection state and the contact surface discharge coefficient, updating the pore water pressure, the shear strength and the strength degradation index; combining the rear edge displacement to calculate the phase error advance time, and according to the phase error advance time, deducing the landslide instability probability, the sliding influence range and the disaster-bearing body dynamic risk value in advance. The present application can identify the hidden degradation of the contact surface when the slope surface displacement has not yet responded significantly, and improve the timeliness and accuracy of the loess-red bed landslide risk deduction under rainfall conditions.

Description

Technical Field

[0001] This invention relates to the field of computer-aided simulation technology, specifically to a method for dynamic simulation of landslide risk in loess red beds based on digital twins. Background Technology

[0002] Loess-red bed landslides typically develop in slope structures where the overlying loess layer is in contact with the underlying red bed, mudstone layer, or weak interlayer. Their stability is influenced by the topographic slope, the degree of development of loess porosity and fissures, the water-blocking effect of the top surface of the red bed, the shear strength of the loess-red bed contact surface, and the rainfall infiltration process. Under rainfall, rainwater can easily be rapidly introduced into the slope body along the cracks at the top of the slope, the tensile cracks at the rear edge, or the unloading cracks, forming local stagnant water near the loess-red bed contact surface. This increases the pore water pressure at the contact surface and reduces the effective normal stress, while also causing a decrease in cohesion and internal friction angle. Consequently, the landslide body evolves from a basically stable state to an unstable or even instable state.

[0003] Current risk assessments for loess-red bed landslides primarily rely on rainfall, stability coefficients, slope displacement, instability probability, or the static exposure relationship of the affected body. While these methods can reflect the hazard and risk level under specific conditions to some extent, they fail to adequately represent the hidden process by which rainfall, after being diverted through rear-edge fissures, forms a water-stagnant wedge at the loess-red bed interface, causing the shear strength of the interface to decay prior to slope displacement. Furthermore, due to the temporal discontinuity between the internal strength degradation of the interface and the response to rear-edge surface displacement, updating the landslide instability probability and sliding impact range only after a significant increase in slope displacement can easily lead to a lag in risk projection. This prevents the spatial, temporal, and dynamic risk values ​​of the affected body from promptly reflecting the expansion and changes in the landslide impact range during rainfall, thus affecting the accuracy of landslide early warning and risk management. Summary of the Invention

[0004] The purpose of this invention is to address the problem in the prior art that the internal strength degradation of the loess-red bed contact surface under rainfall is difficult to characterize in a timely manner, resulting in a lag in the dynamic estimation of landslide instability probability, sliding impact range, and risk of the disaster-bearing body. The invention proposes a dynamic estimation method for loess-red bed landslide risk based on digital twins.

[0005] The technical solution of this invention: a method for dynamic estimation of loess red bed landslide risk based on digital twins, comprising: S1. Obtain topographic profile, loess layer boundary, red layer top surface boundary, loess-red layer contact surface location, rear edge crack location, potential sliding boundary and disaster-bearing body spatial distribution data of the target landslide, construct a digital twin contact surface model of the loess-red layer landslide, and obtain basic data of the contact surface unit; S2. Input the rainfall intensity, and calculate the water retention height of the contact surface unit by combining the trailing edge crack connectivity state and the equivalent discharge coefficient of the contact surface unit in the basic data of the contact surface unit. S3. Update the pore water pressure and shear strength of the contact surface unit according to the water retention height, and calculate the strength degradation index; S4. Based on the strength degradation index and the trailing edge surface displacement, calculate the phase advance time of the strength degradation relative to the trailing edge displacement response. S5. Calculate the landslide instability probability based on the shear strength, and map the landslide instability probability to the predicted instability probability based on the phase shift lead time. S6. Generate the slippage impact range based on the instability probability and the lead time of the phase shift, and calculate the dynamic risk value of the disaster-bearing body by combining the spatial distribution data of the disaster-bearing body.

[0006] Preferably, the basic data of the contact surface unit includes the contact surface unit length, contact surface inclination angle, overburden weight per unit width, initial cohesion, initial internal friction angle, initial pore water pressure, contact surface equivalent discharge coefficient, and trailing edge fracture connectivity status; the trailing edge fracture connectivity status includes a connected state and a disconnected state, wherein a connected state indicates that the corresponding contact surface unit and the trailing edge fracture have a hydraulic connection, and a disconnected state indicates that the corresponding contact surface unit and the trailing edge fracture do not have a hydraulic connection.

[0007] Preferably, calculating the water retention height of the contact surface unit includes: The equivalent infiltration flux into the contact surface unit is determined based on the rainfall intensity and the connectivity of the trailing fractures. The discharge flux of the contact surface unit is determined based on the equivalent discharge coefficient of the contact surface, the water height at the previous moment, and the length of the contact surface unit. Based on the previous water level height, equivalent infiltration flux, discharge flux, and calculation time step, the current water level height is obtained according to the water conservation relationship.

[0008] Preferably, updating the pore water pressure and shear strength of the contact surface unit includes: The pore water pressure at the current moment is determined based on the initial pore water pressure, the specific weight of water, and the current water level. Based on the initial cohesion, initial internal friction angle, and water content degradation relationship, determine the cohesion and internal friction angle at the current moment; The total normal stress is determined based on the component of the overburden weight per unit width in the normal direction of the contact surface and the length of the contact surface element. The shear strength at the current moment is determined according to the Mohr-Coulomb shear strength criterion based on the cohesion, total normal stress, pore water pressure, and internal friction angle at the current moment.

[0009] Preferably, the strength degradation index is calculated, including: The shear strength before the start of rainfall is taken as the initial shear strength, and the ratio of the difference between the initial shear strength and the current shear strength to the initial shear strength is determined as the strength deterioration index.

[0010] Preferably, calculating the phase misalignment advance time includes: For contact surface elements in a connected state, the strength degradation index is weighted and averaged according to the length of the contact surface element to obtain the strength degradation amount of the trailing edge control area. The trailing edge displacement velocity is determined based on the difference in trailing edge surface displacement at adjacent calculation times and the calculation time step. Time-delay correlation calculations were performed on the intensity degradation of the trailing edge control zone and the trailing edge displacement velocity. The candidate advance time corresponding to the maximum correlation was determined as the phase misalignment advance time.

[0011] Preferably, calculating the probability of instability in advance includes: Based on the current shear strength and the length of the contact surface element, determine the resultant force of the anti-slip force along the contact surface according to the unit width section; The resultant force of the sliding force along the contact surface is determined based on the weight of the overlying soil per unit width and the inclination angle of the contact surface. The stability coefficient is determined based on the ratio of the resultant anti-slip force to the resultant sliding force. The probability of landslides under rainfall conditions is determined based on the correlation between the stability coefficient and the probability of landslides. The probability of occurrence of rainfall-induced factors is determined based on the frequency of occurrence of the current rainfall process or the set working conditions, and the probability of landslide instability is determined based on the probability of occurrence of rainfall-induced factors and the probability of landslide under rainfall conditions. By mapping the landslide instability probability to the projected time after the phase shift, the projected instability probability is obtained.

[0012] Preferably, the generated slip influence range includes: The instability probability and the advance time of phase misalignment are input into the digital twin slip range model verified by historical numerical simulation results. The digital twin slip range model uses the maximum slip distance and the landslide instability impact range under rainfall conditions as verification objects to establish the correspondence between the instability probability and the slip impact range, and generates the slip impact range corresponding to the instability time after the advance time of phase misalignment.

[0013] Preferably, the calculation of the dynamic risk value of the disaster-bearing body includes: By overlaying the impact range of the slippage with the spatial distribution data of the disaster-bearing bodies, the spatial probability of disaster-bearing bodies can be determined. Determine the disaster-bearing body's temporal probability of disaster, disaster-bearing body's loss rate, and disaster-bearing body's value based on the disaster-bearing body's category; The dynamic risk value of a disaster-bearing body is determined by multiplying the instability probability, spatial disaster-bearing probability, temporal disaster-bearing probability, loss rate, and value of the disaster-bearing body in advance.

[0014] Preferably, the disaster-bearing body category includes personnel, buildings, and roads; when the disaster-bearing body category is personnel, the disaster-bearing body value is expressed in terms of the number of personnel, and the dynamic risk value of the disaster-bearing body is expressed in terms of the number of people; when the disaster-bearing body category is buildings or roads, the disaster-bearing body value is expressed in terms of economic value, and the dynamic risk value of the disaster-bearing body is expressed in terms of monetary amount.

[0015] Compared with the prior art, the above-mentioned technical solution of the present invention has the following beneficial technical effects: This invention constructs a digital twin contact surface model of loess-red bed landslides, unifying the topographic profile of the target landslide, the loess-red bed contact surface, the rear edge fissures, the potential sliding boundary, and the spatial distribution data of the disaster-bearing body into the same dynamic calculation framework. Furthermore, the loess-red bed contact surface is discretized into multiple contact surface units, enabling continuous expression of rainfall infiltration, fissure diversion, contact surface water retention, pore water pressure increase, and shear strength decay at specific spatial locations. Instead of relying solely on rainfall, slope displacement, or overall stability coefficient for risk assessment, it directly characterizes the hidden degradation process inside the loess-red bed contact surface under rainfall through the progressive updates of water retention height, pore water pressure, cohesion, internal friction angle, and shear strength, thereby improving the physical relevance and dynamic response capability of loess-red bed landslide risk estimation. This invention further obtains the phase advance time of strength degradation relative to the rear displacement response by calculating the time lag correlation between the strength degradation index and the rear displacement velocity. This phase advance time is then used to map the current landslide instability probability to the future projection time. Combined with a verified digital twin slip range model, the slip impact range is generated. Finally, the spatial range, temporal exposure probability, loss rate, and value of the disaster-bearing body are superimposed to obtain the dynamic risk value of the disaster-bearing body. This invention can identify the instability risk caused by the premature attenuation of the loess-red bed contact surface strength before the slope surface displacement increases significantly. This allows the landslide instability probability, slip impact range, and disaster-bearing body risk value to be updated synchronously with the rainfall process, reducing the risk projection lag and improving the timeliness and accuracy of landslide early warning, personnel evacuation, road control, and building risk prevention under rainfall conditions. Attached Figure Description

[0016] Figure 1 This is a flowchart of the dynamic simulation method for loess red layer landslide risk based on digital twin proposed in this invention; Figure 2 This is a schematic diagram of the water-stagnant wedge and strength degradation at the contact surface under the trailing edge crack flow guidance proposed in this invention. Detailed Implementation

[0017] Example 1, as Figure 1 As shown in the figure, this embodiment provides a method for dynamic estimation of the risk of loess red layer landslide based on digital twins. It is used to continuously calculate and dynamically express the process of loess red layer landslide under rainfall, including back edge fissure diversion, water retention at the loess-red layer contact surface, deterioration of the shear strength of the contact surface, hysteresis response of slope surface displacement, early estimation of landslide instability probability, and dynamic risk update of the disaster-bearing body. To avoid confusion between the calculation time step, candidate advance time, and phase shift advance time, in this embodiment, the calculation time step is uniformly denoted as... The unit is seconds (s); the candidate advance time is uniformly denoted as... The unit is seconds (s); the phase misalignment advance time is uniformly denoted as... The unit is seconds; the maximum candidate advance time is uniformly denoted as... The unit is seconds (s); the length of the sliding calculation window is uniformly denoted as... The unit is seconds (s); the advance simulation time will be uniformly represented as... In this embodiment, all time intervals involving recursive calculations of adjacent moments are calculated using... All time parameters involving the lead time of strength degradation relative to the trailing edge displacement response are adopted. , or This makes the physical meaning of each time parameter independent of each other.

[0018] During implementation, the area where the target landslide is located is used as the calculation object. Topographic profiles, loess layer boundaries, red bed top surface boundaries, loess-red bed contact surface locations, rear edge fissure locations, and spatial distribution data of the disaster-bearing bodies are obtained through UAV oblique photography, lidar measurement, total station measurement, engineering geological mapping, borehole exposure, trench exposure, groundwater monitoring, rear edge fissure investigation, and existing geological survey data. All of the above data are then unified under the same plane coordinate system and elevation datum to form three-dimensional topographic data, stratigraphic interface data, fissure line data, potential slip boundary data, and disaster-bearing body layer data of the target landslide. Among these, the spatial distribution data of the disaster-bearing bodies includes locations or areas of human activity, building outlines, road centerlines, road widths, building values, road values, and the number of people. A digital twin contact surface model of a loess-red layer landslide is constructed in a computer server or edge computing terminal. In this embodiment, the digital twin contact surface model of the loess-red layer landslide refers to a virtual calculation model formed in a computer based on the topographic profile of the target landslide, the boundary of the loess layer, the top boundary of the red layer, the location of the loess-red layer contact surface, the location of the rear edge fissures, the potential sliding boundary, and the spatial distribution data of the disaster-bearing body. This model can be updated synchronously with rainfall input, the water retention state of the contact surface, the shear strength of the contact surface, the rear edge surface displacement, and the risk status of the disaster-bearing body. This model is not only used for three-dimensional display, but also for calculating the water retention height, updating the pore water pressure, updating the shear strength, calculating the phase difference advance time, calculating the instability probability in advance, and calculating the dynamic risk value of the disaster-bearing body within each calculation time step. The loess-red bed contact surface is divided into multiple contact surface units along the rear-to-front direction, denoted as the i-th contact surface unit. In this embodiment, the contact surface unit refers to the calculation segment formed by discretizing the loess-red bed contact surface according to the main sliding direction of the landslide. Each contact surface unit corresponds to a specific spatial location and an independent mechanical calculation object. Each contact surface unit includes at least the contact surface unit length, contact surface dip angle, overburden weight per unit width, initial cohesion, initial internal friction angle, initial pore water pressure, contact surface equivalent discharge coefficient, and rear edge fracture connectivity. By setting contact surface units, the local water retention and local strength deterioration caused by rainfall infiltration can be mapped to the specific location of the loess-red bed contact surface, avoiding the local deterioration process being expressed only by the overall landslide stability coefficient. Each contact surface unit All correspond to the length of the contact surface unit. Contact surface inclination angle Weight of overlying soil per unit width Initial cohesion Initial internal friction angle Initial pore water pressure Equivalent discharge coefficient of contact surface Connected state with trailing edge fissure ; Length of contact surface unit The contact angle is determined by the arc length (in meters) of the broken line or curve of the loess-red bed contact surface within the contact surface unit; The weight of the overlying soil per unit width is calculated from the elevation difference between the two endpoints of the contact surface unit and the horizontal projection distance. The initial cohesion was calculated from the loess body area per unit width profile above the contact surface unit, the natural unit weight of the loess, and the horizontal projection range corresponding to the contact surface unit, with units of kN / m. and initial internal friction angle The initial pore water pressure is determined by a combination of direct shear tests, circumferential shear tests, field inversion results, or parameters from adjacent engineering surveys of disturbed or undisturbed loess-red bed contact samples. Determined by groundwater level monitoring values, water pressure gauge monitoring values, or steady-state seepage inversion results before rainfall begins; equivalent discharge coefficient at the contact surface. Determined by permeability tests near the loess-red bed contact surface, in-situ rainfall infiltration tests, or historical rainfall displacement inversion; In this embodiment, the equivalent discharge coefficient of the contact surface refers to the equivalent parameter obtained by comprehensively converting the drainage capacity near the loess-red bed contact surface, the water-blocking effect of the top surface of the red bed, the drainage capacity along the contact surface, and the influence of local seepage channels. Its unit is m / s. The equivalent discharge coefficient of the contact surface is not the same as the permeability coefficient of a single rock and soil body. Instead, it is a parameter used to calculate the dissipation rate of the stagnant water height of the contact surface unit. The larger this parameter is, the easier it is for the stagnant water formed by the contact surface unit to be discharged along the contact surface or the adjacent seepage channel.

[0019] trailing edge gap connectivity This is used to indicate whether there is hydraulic connectivity between the i-th contact surface unit and the trailing edge fracture. In this embodiment, the trailing edge fracture connectivity state is a binary state quantity used to describe whether rainwater can enter the corresponding contact surface unit through the trailing edge fracture, fracture zone, unloading zone, or weak interlayer. The trailing edge fracture connectivity state only characterizes the hydraulic connectivity relationship and does not characterize whether the contact surface unit itself has become unstable. When the trailing edge fracture connectivity state is connected, the contact surface unit participates in the calculation of the stagnant water height under the guidance of the trailing edge fracture. When the trailing edge fracture connectivity state is not connected, the contact surface unit does not directly receive the flow replenishment from the trailing edge fracture. When the trailing edge fracture extends spatially to the contact surface unit, or forms a continuous water-conducting channel with the contact surface unit through the fracture zone, unloading zone, or weak interlayer, it is determined that... When there is no continuous water-conducting channel between the trailing edge fracture and the contact surface unit, or when the middle is blocked by a low-permeability layer, it is determined that... ; This yields the basic data of the contact surface element, which contains... , and In subsequent calculations of the water retention height of the contact surface element, the basic data of the contact surface element... , , , , and The spatial distribution data of the disaster-bearing body is subsequently used to update pore water pressure, shear strength, and landslide instability probability. This data is then used to calculate the dynamic risk value of the disaster-bearing body, thus establishing a correspondence between the modeling data, mechanical calculation data, and risk calculation data.

[0020] During rainfall, rainfall intensity is input in real time through rain gauges, weather stations, or radar rainfall grids. The rainfall intensity was uniformly converted to m / s; the equivalent infiltration flux into each contact surface unit was determined according to the trailing edge fracture connectivity. The calculation method is as follows: ; in, Let be the equivalent infiltration flux of the i-th contact surface element at time t, in m / s; The trailing edge fissure is connected and is dimensionless. Let be the rainfall intensity at time t, expressed in m / s. This calculation method is used because connected contact surface elements can receive rapid rainfall infiltration from trailing edge fissures, while disconnected contact surface elements do not directly receive rainfall from trailing edge fissures. As a switching quantity for rainfall entering the loess-red bed contact surface, it can express the control effect of the back-edge fissure flow on the formation of water stagnation at the contact surface; Based on the equivalent discharge coefficient of the contact surface The water level at the previous moment and contact surface unit length Determine the discharge flux of the contact surface unit. The calculation method is as follows: ; in, Let f be the flux discharged along the contact surface by the i-th contact surface element at time t, in m / s; The equivalent discharge coefficient of the contact surface is expressed in m / s. The water level of the contact surface element at the previous calculation time is in meters. The length of the contact surface unit is given in meters (m). This calculation method is used because a greater stagnant water height results in a larger hydraulic gradient along the loess-red bed contact surface, leading to a greater discharge flux. Conversely, a larger contact surface unit length results in a smaller hydraulic gradient per unit length. Therefore, the length of the contact surface unit is given as... This indicates that the drainage driving gradient can reflect the gradual discharge process of water stuck at the contact surface; Then, based on the water level at the previous moment... Equivalent infiltration flux Discharge flux and calculation time step The current water level is obtained according to the law of conservation of water volume. The calculation method is as follows: ; in, Let be the water level height of the i-th contact surface element at time t, in meters. The time step is calculated in seconds. The reason for using this calculation method is that the current water level is determined by the water level at the previous moment, the current rainfall inflow, and the current discharge. This method can continuously reflect the accumulation and dissipation of water wedges formed after rainfall enters the loess-red bed contact surface through the rear edge fissures. In this embodiment, a water-retention wedge refers to a localized water-retention area formed above the loess-red bed contact surface after rainfall is introduced into the loess layer through rear-edge fissures. This area is formed due to the relative water resistance of the red bed top surface and the limited drainage capacity of the contact surface. The water-retention wedge is not required to have a strictly wedge-shaped geometry; rather, it characterizes the hydraulic state where the water retention height above the contact surface is differentially distributed along the rear-to-front direction. The water retention height is used to quantify the local water head size of the water-retention wedge at each contact surface unit and serves as a direct input for updating pore water pressure. Figure 2 As shown, rainwater infiltrates downwards along cracks at the rear edge of the slope and flows into the loess layer. Near the loess-red bed contact surface, it is constrained by the relative water-blocking effect of the red bed and the limited drainage capacity of the contact surface, forming a water-stagnant wedge. The water-stagnant wedge increases the water-stagnant height of the corresponding contact surface unit, which in turn causes an increase in pore water pressure, a decrease in effective normal stress, and a decrease in the shear strength of the loess-red bed contact surface. Since the internal strength deterioration of the contact surface occurs before the significant response of the slope surface displacement, the surface displacement at the rear edge exhibits a delayed response. By identifying this phase relationship of strength decreasing first and displacement appearing later, the probability of landslide instability and the range of landslide influence can be predicted in advance.

[0021] According to the height of the water level Update the pore water pressure of the i-th contact surface element. The calculation method is as follows: ; in, Let be the pore water pressure of the i-th contact surface element at time t, in kPa; This represents the initial pore water pressure, expressed in kPa. The specific gravity of water, measured in units of 1000 kJ / m². ; The water level is the stagnant water height, in meters (m). Equal to 1 kPa, therefore This corresponds to the additional pore water pressure caused by the increase in the height of the water retention. The reason for using this calculation method is that the rainwater introduced by the rear edge cracks will directly increase the pore water pressure of the contact surface after forming water retention above the loess-red layer contact surface. The increase in pore water pressure will reduce the effective normal stress, which is the direct mechanical cause of the deterioration of the shear strength of the contact surface. Based on initial cohesion Initial internal friction angle And the relationship between water content degradation and cohesion at the current moment and the friction angle at the current moment In this embodiment, the water content degradation relationship refers to the relationship between the cohesion and internal friction angle of the loess-red bed contact surface material under different water content conditions. This relationship is used to convert the water softening effect caused by the water retention height into the changes in cohesion and internal friction angle at the current moment. The water content degradation relationship can be obtained through direct shear test, ring shear test, triaxial test, field inversion or historical rainfall conditions under different water content conditions. Its function is to avoid characterizing the strength decay only by the change in pore water pressure, thereby simultaneously reflecting the decrease in strength parameters caused by the water softening of the contact surface material itself. As one specific implementation method, the ratio of the water retention height to the thickness of the contact surface is determined as the degree of water content degradation. The calculation method is as follows: ; in, The degree of water content degradation of the i-th contact surface unit; The contact surface influence thickness of the contact surface unit is in meters. It can be determined by the thickness of the weak interlayer revealed by the borehole, the thickness of the weathered zone on the top surface of the red bed, the thickness revealed by the field trench, or the thickness of the contact surface calculated in the numerical model. The reason for using this calculation method is that the closer the water retention height is to the contact surface influence thickness, the more fully the loess-red bed contact surface is affected by water softening. When the water retention height reaches or exceeds the contact surface influence thickness, the degree of deterioration is treated as fully deteriorated. Current moment cohesion and the friction angle at the current moment Determined according to the following formula: ; ; in, Let be the cohesive force of the i-th contact surface element at time t, in kPa; The residual cohesion of the i-th contact surface unit is expressed in kPa. The initial cohesion of the i-th contact surface unit is expressed in kPa. The water content degradation coefficient represents the cohesion. Let be the internal friction angle of the i-th contact surface element at time t, in degrees; The residual internal friction angle of the i-th contact surface element is expressed in degrees. Let be the initial internal friction angle of the i-th contact surface element, in degrees; The internal friction angle is the water content deterioration coefficient, in degrees; , , and All of these can be determined by indoor direct shear tests, ring shear tests, field inversion, or calibration of historical landslide conditions under different water content states. The reason for adopting the above-mentioned water content deterioration relationship is that after the loess-red layer contact surface is wetted by water, both the cohesion and the internal friction angle will decrease with the increase of water content. However, the decrease is constrained by the residual strength. Using the lower limit of the residual strength can avoid calculating negative strength or excessively low strength that does not conform to the meaning of geotechnical mechanics.

[0022] Based on the weight of the overlying soil per unit width The component in the normal direction of the contact surface and the length of the contact surface element. Determine the total normal stress The calculation method under the condition of unit width section is as follows: ; in, Let be the total normal stress of the i-th contact surface element at time t, in kPa; The weight of the overlying soil per unit width of the i-th contact surface element is expressed in kN / m. The contact surface inclination angle; The length of the contact surface element is given in meters (m). This calculation method is used because... The component of the weight of the overlying soil in the normal direction of the contact surface is given by dividing it by the length of the contact surface element to obtain the total normal stress acting on the contact surface under the condition of unit width. Then, based on the cohesion at the current moment Total normal stress Current pore water pressure and the friction angle at the current moment Determine the shear strength at the current moment according to the Mohr-Coulomb shear strength criterion. The calculation method is as follows: ; in, Let be the shear strength of the i-th contact surface element at time t, in kPa; The cohesive force at the current moment, expressed in kPa; The total normal stress is expressed in kPa. The pore water pressure at the current moment is expressed in kPa. The internal friction angle is the value at the current moment. This calculation method is used because the Mohr-Coulomb shear strength criterion can simultaneously reflect the contributions of cohesion, effective normal stress, and internal friction angle to the shear capacity of the contact surface. The water retention height reduces the effective normal stress through pore water pressure and decreases cohesion and internal friction angle through the water content degradation relationship, thus expressing the process of the initial attenuation of the strength of the loess-red bed contact surface under rainfall conditions. The reason is that when the interfacial water pressure at the contact surface exceeds the total normal stress, the effective compressive stress at the contact surface no longer participates in the friction resistance calculation as a negative value, thus avoiding obtaining a negative friction resistance that does not conform to physical meaning.

[0023] Using the shear strength before the onset of rainfall as the initial shear strength, the ratio of the difference between the initial shear strength and the current shear strength to the initial shear strength is determined as the strength degradation index. The calculation method is as follows: ; in, Let be the strength degradation index of the i-th contact surface element at time t; This refers to the initial moment before the rainfall begins. Let be the shear strength of the i-th contact surface element at the initial moment, in kPa; Let be the shear strength of the i-th contact surface element at time t, in kPa. The reason for using this calculation method is that the strength degradation index can convert the reduction in shear strength of different contact surface elements into a dimensionless proportion, which is convenient for subsequent unified comparison and weighted summarization of contact surface elements with different lengths and different physical parameters. In this embodiment, the strength degradation index is a dimensionless index used to characterize the degree of decrease in the current shear strength of the contact surface unit relative to the initial shear strength before the start of rainfall. The strength degradation index only reflects the relative degree of shear strength decay and is not directly equivalent to the probability of landslide instability. Subsequently, it participates in the early simulation through the strength degradation amount of the trailing control zone and the phase shift advance time to identify whether the strength decay inside the contact surface occurs before the trailing displacement response.

[0024] To calculate the phase advance time of strength degradation relative to the trailing edge displacement response, this embodiment performs a weighted average of the strength degradation index on the contact surface elements in the connected state according to the contact surface element length, thereby obtaining the strength degradation amount in the trailing edge control region. The calculation method is as follows: ; in, The strength degradation amount in the trailing edge control zone; n is the total number of contact surface elements; The trailing edge crack is in a connected state; Let be the strength degradation index of the i-th contact surface element; The length of the i-th contact surface element is in meters. The reason for using this calculation method is that only the contact surface elements that are in hydraulic communication with the trailing edge fracture can directly reflect the hidden strength degradation caused by the flow of the trailing edge fracture. Moreover, the larger the length of the contact surface element, the greater its contribution to the overall strength state of the trailing edge control area. Therefore, using length-weighted average can avoid the excessive influence of the outlier values ​​of short elements on the overall degradation judgment. In this embodiment, the trailing edge control zone refers to the contact surface area jointly formed by contact surface units in the trailing edge crack connectivity state; the strength degradation of the trailing edge control zone is a comprehensive degradation index obtained by weighting the strength degradation index of each contact surface unit in the trailing edge control zone according to the length of the contact surface unit, which is used to characterize the overall strength attenuation within the trailing edge crack flow control area; this index is used for subsequent time-delay correlation calculation with the trailing edge displacement velocity to identify whether the strength degradation precedes the slope surface displacement response; when At that time, the intensity degradation of the trailing edge control zone is calculated according to the above length-weighted method. ;when When this occurs, it indicates that there are no contact surface elements in the current digital twin contact surface model that form a hydraulic connection with the trailing edge fracture. In this case, Set it to 0, and advance the subsequent phase shift time. Setting it to 0 indicates that the pre-deterioration process of the loess-red bed contact surface strength caused by the backward fissure conduction was not identified at the current calculation time. The reason for adopting this boundary treatment method is that the strength deterioration of the backward control zone is used to characterize the attenuation of the contact surface strength in the backward fissure conduction control zone. When there are no connected contact surface elements, there is no physical basis for calculating the length-weighted deterioration. Therefore, setting 0 as the calculation result without backward fissure conduction deterioration can avoid the model interruption caused by the denominator being zero.

[0025] Simultaneous acquisition of trailing edge surface displacement The rear edge surface displacement can be obtained from global navigation satellite system monitoring points, crack gauges, surface displacement gauges, synthetic aperture radar interferometry results, or UAV image matching results, and uniformly projected onto the main sliding direction of the landslide; the rear edge displacement velocity is determined based on the difference in rear edge surface displacement at adjacent calculation times and the calculation time step. The calculation method is as follows: ; in, The trailing edge displacement velocity is expressed in m / s. The displacement of the trailing edge surface at time t is expressed in meters. This represents the rear edge surface displacement at the previous calculation time, in meters. The time step is calculated in seconds. This calculation method is used because displacement velocity is better than cumulative displacement in reflecting the response characteristics of the landslide rear edge as it transforms from a stable state to an accelerated deformation state in a short period of time. In length Within the sliding calculation window, the intensity degradation of the trailing edge control region is... and trailing edge displacement velocity Perform time-delay related calculations; the length of the sliding calculation window. The determination is made as follows: When there is a historical rainfall monitoring process for the target landslide, the shortest continuous time period that can cover the transition from a steady change to an accelerated change in the rear edge displacement velocity is selected as the [missing information]. When historical rainfall monitoring data is lacking for the target landslide, Take at least three computation time steps And not less than the shortest time required for trailing edge displacement monitoring data to complete one effective change identification; the maximum candidate advance time The determination is made as follows: When historical rainfall events or numerical simulation samples exist, the time difference between the peak time of intensity degradation in the trailing control zone and the peak time of the trailing displacement velocity is extracted, and the largest positive time difference is taken as the threshold. When no historical rainfall events or numerical simulation samples exist, one-third of the current rainfall duration is taken as the baseline. Candidate advance time Starting from 0, according to the calculation time step Incrementing values ​​until a value is reached. ; For each candidate advance time Calculate the strength degradation of the trailing edge control zone. The trailing edge displacement velocity after translation according to the candidate advance time Correlation coefficient between them: ; in, advance time for candidates The corresponding correlation coefficient; For sliding calculation windows; This represents the average intensity degradation of the trailing edge control zone within the sliding calculation window. The mean trailing edge displacement velocity after translation according to the candidate advance time within the sliding calculation window; when any standard deviation term in the correlation coefficient calculation formula is 0, the correlation coefficient corresponding to the candidate advance time is set to 0, so as to avoid the correlation coefficient being unable to be calculated due to the lack of data fluctuation; The candidate advance time corresponding to the maximum correlation coefficient is determined as the phase misalignment advance time:

[0026] in, The phase advance time of strength degradation relative to the trailing edge displacement response, in seconds; The candidate lead time is expressed in seconds (s). In this embodiment, the phase advance time refers to the advance in time between the strength degradation of the control zone and the displacement velocity response at the rear edge. The phase advance time is used to characterize the time difference between the occurrence of strength degradation inside the loess-red layer contact surface and the lack of synchronous and significant response of the surface displacement at the rear edge. When the phase advance time is 0, it indicates that no reliable relationship between strength degradation and displacement is identified in the current calculation window. When the phase advance time is greater than 0, it indicates that the strength degradation inside the contact surface has a calculable advance relative to the displacement response at the rear edge, and serves as the time mapping basis for subsequent advance estimation of instability probability and slippage influence range. The lower limit of correlation is determined as follows: When there is a set of historical rainfall events without instability, the maximum correlation coefficient corresponding to each candidate advance time in that historical rainfall event is calculated, and this maximum value is used as the lower limit of correlation; when there are multiple sets of historical rainfall events without instability and the number of sets is not less than 5, the maximum correlation coefficient corresponding to each candidate advance time in each set of historical rainfall events without instability is calculated, and the 95th percentile of all maximum values ​​is used as the lower limit of correlation; when there are multiple sets of historical rainfall events without instability but the number of sets is less than 5, the maximum value among all maximum values ​​is used as the lower limit of correlation; when there is no historical rainfall event without instability, 0.6 is used as the lower limit of correlation; if If it falls below the lower limit of correlation, then... Setting it to 0 indicates that no reliable phase degradation relationship of initial strength decrease followed by subsequent displacement has been identified within the current calculation window; the reason for using this method is that the sliding calculation window... The length of data used for phase error detection and the maximum candidate advance time are determined. The searchable time range for intensity degradation preceding trailing edge displacement response is limited, and the lower correlation limit is used to exclude spurious misphase relationships caused by random fluctuations.

[0027] Based on the current shear strength and contact surface unit length Determine the resultant anti-slip force along the contact surface based on the unit width profile. The calculation method is as follows: ; in, The resultant force of the anti-sliding force along the loess-red bed contact surface at time t is expressed in kN / m. Let be the shear strength of the i-th contact surface element at time t, in kPa; The length of the i-th contact surface element is in meters. The reason for using this calculation method is that, under the condition of a unit width section, the shear strength of the contact surface multiplied by the length of the contact surface element can be used to obtain the anti-slip force that the element can provide. The anti-slip forces of each contact surface element can be added together to obtain the overall anti-slip capacity along the contact surface. Based on the weight of the overlying soil per unit width and contact surface inclination angle Determine the resultant force of the sliding force along the contact surface The calculation method is as follows: ; in, The resultant force of the sliding force along the loess-red bed contact surface at time t is expressed in kN / m. The weight of the overlying soil per unit width of the i-th contact surface element is expressed in kN / m. Let be the contact surface inclination angle of the i-th contact surface element; the reason for using this calculation method is that... The sliding component generated by the weight of the overlying soil in the i-th contact surface unit along the contact surface direction is represented by the sum of the sliding components of each unit, which can characterize the driving force for the landslide to slide along the loess-red layer contact surface. Based on the resultant force of anti-slip force Combined with the downward force The ratio determines the stability coefficient. The calculation method is as follows: ; in, Let t be the stability coefficient of the target landslide at time t. The reason for using this calculation method is that the stability coefficient can characterize the relative magnitude of the anti-sliding effect to the sliding effect. When the resultant force of the anti-sliding force decreases relative to the resultant force of the sliding force, the possibility of the landslide entering an unstable state increases. Based on stability coefficient The correspondence between rainfall and landslide probability determines the probability of landslide occurrence under rainfall conditions. As a specific implementation method, a stability coefficient-slip probability grading table is used to determine the slip probability under rainfall conditions. This stability coefficient-slip probability grading table can be established based on landslide stability state classification standards, historical landslide samples, local landslide risk assessment specifications, or numerical simulation samples. For example, a stability coefficient less than 1.00 corresponds to a slip probability range of 0.50 to 1.00; a stability coefficient greater than or equal to 1.00 and less than 1.10 corresponds to a slip probability range of 0.40 to 0.50; a stability coefficient greater than or equal to 1.10 and less than 1.20 corresponds to a slip probability range of 0.30 to 0.40; and a stability coefficient greater than or equal to 1.20 and less than or equal to 2.00 corresponds to a slip probability range of 0.10 to 0.30. When the value falls within a certain stability coefficient range, linear interpolation is used to determine the stability coefficient within that range. When the stability coefficient is lower than the minimum boundary of the classification table, the highest slip probability of the classification table is taken; when the stability coefficient is higher than the maximum boundary of the classification table, the lowest slip probability of the classification table is taken. The reason for adopting this method is that the relationship between the stability coefficient and the slip probability is not a single-point jump relationship, but can be mapped between intervals according to the stable state classification. Linear interpolation can form continuous probability results between the boundaries of the classification table and reduce the uncertainty caused by the introduction of additional calibration parameters.

[0028] Determine the probability of occurrence of rainfall-inducing factors based on the frequency of current rainfall events or the set operating conditions. When using historical rainfall frequency, it can be obtained from statistics of current cumulative rainfall, duration, and return period; when using early warning setting conditions, it can be determined according to the setting probability of that rainfall condition in the early warning system; and then based on the occurrence probability of rainfall-inducing factors. Probability of landslides under rainfall conditions Determine the probability of landslide instability based on the lead time of phase shift. By mapping the landslide instability probability to the predicted time, we obtain the predicted instability probability. The calculation method is as follows: ; in, To advance the simulation The probability of landslide instability at any given moment; Let t be the probability of the rainfall-inducing factor occurring at time t; The probability of slippery surfaces under rainfall conditions; The phase advance time is expressed in seconds. This calculation method is used because the probability of rainfall-induced factors represents the likelihood of the current rainfall condition occurring, while the probability of landslide instability under rainfall conditions represents the likelihood of the landslide becoming unstable due to a decrease in the stability coefficient. Multiplying the two yields the probability of rainfall-induced landslide instability, which is then mapped to... It can always demonstrate the early prediction effect of the contact surface strength degradation preceding the significant response of the trailing edge displacement; In this embodiment, the predicted instability probability refers to the landslide instability probability calculated based on the shear strength of the contact surface at the current moment, and then mapped to the predicted future moment based on the phase shift advance time. The predicted instability probability is not an additional independent probability, but rather a forward mapping of the instability trend already revealed by the current contact surface strength degradation to... The results at any given moment are used to update the slip impact range and dynamic risk value of the disaster-bearing body in advance, before the trailing edge displacement has fully responded. Based on the predicted probability of instability and phase advance time The digital twin slip range model, verified by historical numerical simulation results, is used to generate the predicted timing. The corresponding slippage impact range; In this embodiment, the digital twin slippage range model refers to the calculation sub-model in the digital twin contact surface model of the loess red layer landslide, which is used to convert the predicted instability probability and the advance time of phase misalignment into the maximum slippage distance and slippage impact range; This model is not an unconstrained black box model. Its parameters are verified by historical landslide inversion results, rainfall condition numerical simulation results, or field monitoring playback results. Its output slippage impact range is jointly constrained by the maximum allowable slippage distance, potential slippage boundary, main slippage direction of the landslide, topographic slope aspect, valley boundary, and historical numerical simulation verification boundary; The digital twin slip range model uses historical landslide inversion results, rainfall simulation results, or field monitoring playback results as verification samples. Each verification sample includes rainfall conditions, contact surface strength parameters, stability coefficient, instability probability, phase shift lead time, maximum slip distance, and landslide instability impact range. During implementation, the sample instability probability is first extracted from each verification sample. Sample phase misalignment advance time and the maximum slip distance of the sample Where j is the verification sample number; then the sample instability probability is... and sample phase advance time As an input variable, the maximum slip distance of the sample As a target value, the probability impact index Phase advance time correction factor Perform least-squares fitting; where the minimum influence distance in the low-probability slip state is... The minimum slip distance is taken from the low-probability slip verification sample, and the maximum slip distance is obtained from the rainfall condition verification. The maximum slip distance and the maximum candidate lead time are taken from the rainfall condition verification sample. The maximum candidate advance time determined in the phase advance time calculation step is used; When the number of verification samples is insufficient to fit simultaneously and season The maximum slip distance is calculated solely based on the pre-calculated probability of instability; when only one rainfall condition verification sample exists, the maximum slip distance is taken as... The maximum slip distance corresponding to the sample under this rainfall condition was used as the verification value. When new historical numerical simulation results, in-situ slip range inversion results, or in-situ monitoring playback results are added, the data should be re-analyzed. and The reason for adopting the above-mentioned insufficient sample handling rule is that fitting multiple parameters at the same time when the verification sample is insufficient will lead to unstable fitting results. By simplifying the parameters, the sliding range can still be deduced based on the existing verification samples, and the model parameters can be gradually corrected as the subsequent samples are added.

[0029] After completing the parameter verification, the instability probability in advance is mapped to the maximum slip distance corresponding to the in advance time. : ; in, The maximum slip distance corresponding to the predicted time is given in advance, and the unit is meters. The minimum influence distance under low-probability slip conditions, in meters; The maximum slip distance obtained from the rainfall condition verification is expressed in meters. To anticipate the probability of instability; The probability influence index; This is the correction factor for phase misalignment advance time; The phase advance time is expressed in seconds (s). The maximum candidate lead time is expressed in seconds; when or season The entire value is set to 0 to avoid calculations involving division by zero. In this case, the maximum slip distance is dominated by the instability probability calculated in advance. To avoid the maximum slip distance exceeding the physical boundary supported by historical numerical simulations or in-situ inversions, Apply upper limit constraints: ; in, This represents the maximum slip distance after the upper limit constraint, in meters. The maximum allowable slip distance, expressed in meters, is obtained from historical numerical simulation results, in-situ slip range inversion results, or verification under typical rainfall conditions. Subsequent calculations of the slip influence range will use this maximum allowable slip distance. As the maximum sliding distance; The reason for adopting the above verification method is that the higher the probability of instability is predicted in advance, the more likely the landslide impact range is to extend to a greater distance; the more obvious the advance time of phase misalignment, the earlier the internal deterioration of the contact surface occurs relative to the slope surface displacement response, and the impact range of the sliding needs to be expanded in advance when predicting risks.

[0030] Maximum slip distance after obtaining the upper limit constraint Then, using the potential sliding boundary as the starting boundary, the movement is pushed outward along the main sliding direction of the landslide. An initial slip distance envelope is formed; then, based on topographic aspect, valley boundaries, abrupt topographic changes at the slope toe, and historical numerical simulation verification boundaries, the initial slip distance envelope is trimmed and corrected to obtain the predicted time. Corresponding slippage influence range When the initial slip distance envelope exceeds the historical numerical simulation verification boundary, the historical numerical simulation verification boundary is used as the outer boundary; when the initial slip distance envelope does not reach the historical numerical simulation verification boundary, the initial slip distance envelope is used as the outer boundary; when the slip direction is constrained by the valley boundary, the slip influence range expands along the valley direction and does not extend outwards from the watershed; the slip influence range... It is stored in the digital twin contact surface model as a planar layer and serves as the input for subsequent calculation of the disaster-bearing body's spatial disaster probability; In this embodiment, the landslide impact range refers to the spatial area that may be affected by the target landslide, accumulation, or coverage at the time of advance simulation. The landslide impact range is expressed in the form of a planar layer, which includes not only the area where the potential landslide body is located, but also the slope toe area or valley extension area that may be covered after extrapolating along the main sliding direction based on the maximum landslide distance. The landslide impact range is used to superimpose with the spatial range of the disaster-bearing body and directly participates in the calculation of the disaster-bearing body spatial disaster probability. The area affected by the slippage By overlaying the data on the spatial distribution of disaster-bearing bodies, the spatial probability of disaster bearing bodies is determined. To avoid confusion between the spatial range and value of disaster-bearing bodies, in this embodiment, the spatial range of the k-th type of disaster-bearing body is denoted as... The disaster-bearing value of the k-th type of disaster-bearing body is denoted as ;in, Used to characterize areas of human activity, building outlines, road buffer zones, or other disaster-bearing spatial extent. Used to characterize the number of people, the economic value of buildings, the economic value of roads, or the value of road repair; In this embodiment, the spatial range of the disaster-bearing body Used to describe the geometric distribution range of the k-th type of disaster-bearing body in space, such as the area of ​​human activity, building outline, road buffer zone, or radius of human activity points; disaster-bearing body value Used to describe the quantity or economic value of the k-th type of disaster-bearing body; spatial extent of the disaster-bearing body. Only involved in spatial overlay and spatial disaster probability calculation, the value of the disaster-bearing body It is only used in the calculation of dynamic risk values ​​of disaster-bearing bodies, and the two should not be used interchangeably; When the k-th type of disaster-bearing body is a spatial planar disaster-bearing body, the spatial disaster-bearing probability of the disaster-bearing body is determined according to the slip influence range. Spatial range of the disaster-bearing body The overlapping area ratio is determined as follows: ; in, Let be the spatial disaster risk probability of the kth type of disaster-bearing body at the pre-simulated time. To anticipate the extent of the slippage impact at the specified time; The spatial extent of the k-th type of disaster-bearing body; The area is the overlapping portion of the slip-affected area and the spatial extent of the k-th type of disaster-bearing body, in units of... ; The area of ​​the kth type of disaster-bearing body is expressed in units of... ; When the disaster-bearing body is a linear object such as a road, the road centerline is converted into a planar buffer zone according to the road width, and this planar buffer zone is used as... The spatial probability of disaster bearing is then calculated; when the disaster-bearing body is a personnel location, a planar activity area is formed based on the radius of personnel activity, and this planar activity area is used as... The spatial probability of disaster is then calculated; when there is no data on the activity radius of personnel locations, if the personnel locations are located within the slippage influence range... If the location is within the slippage influence area, the probability of spatial disaster is taken as 1. In addition, the spatial probability of disaster is taken as 0. The reason for adopting this calculation method is that the disaster-bearing body is not necessarily affected by the disaster as long as it is located within the administrative area where the landslide is located. Instead, the actual possibility of disaster needs to be determined based on the superposition relationship between the landslide impact range and the spatial location of the disaster-bearing body. Determine the probability of disaster exposure based on the type of disaster-bearing body. Disaster-bearing body loss rate and the value of disaster-bearing bodies The disaster-bearing body categories include personnel, buildings, and roads. When the disaster-bearing body category is personnel, the temporal probability of disaster bearing is determined by the probability that personnel are located in the corresponding area at the pre-simulated time, which can be determined by the resident population, construction worker shift schedules, road traffic flow, or manual inspection records. The disaster-bearing body loss rate can be determined based on the slip intensity, evacuation time, and personnel exposure status. The disaster-bearing body value is expressed in terms of the number of personnel, and the dynamic risk value of the disaster-bearing body is expressed in terms of the number of people. When the disaster-bearing body category is buildings, the temporal probability of disaster bearing can be taken as the probability that the building is in use at that time. The loss rate of a disaster-bearing body can be determined based on the building structure type, the thickness of the slip accumulation, the impact intensity, and the building's damage resistance level. The value of the disaster-bearing body is expressed as the economic value of the building, and the dynamic risk value of the disaster-bearing body is expressed in monetary terms. When the disaster-bearing body category is a road, the temporal probability of disaster-bearing body can be determined based on the traffic flow, traffic period, or road openness at the pre-simulated time. The loss rate of the disaster-bearing body can be determined based on the slip coverage length, coverage thickness, road grade, and repair cost ratio. The value of the disaster-bearing body is expressed as the economic value or repair value of the road, and the dynamic risk value of the disaster-bearing body is expressed in monetary terms.

[0031] The dynamic risk value of the k-th type of disaster-bearing body is determined by multiplying the instability probability, spatial disaster-bearing probability, temporal disaster-bearing probability, loss rate, and value of the disaster-bearing body in advance. ; in, This represents the dynamic risk value of the k-th type of disaster-bearing body at the time of advance simulation; To anticipate the probability of instability; The spatial disaster-bearing probability of the kth type of disaster-bearing body; The probability of disaster bearing capacity over time for the kth type of disaster-bearing body; The loss rate of the kth type of disaster-bearing body; The value of the k-th type of disaster-bearing body; when the k-th type of disaster-bearing body is personnel, Expressed in terms of number of people, Expressed in terms of the number of people; when the k-th type of disaster-bearing body is a building or road, Expressed in monetary terms, Expressed in monetary terms; the reason for adopting this calculation method is that landslide risk is not only determined by whether the landslide is unstable, but also by whether the affected body is within the scope of the landslide, whether it is exposed at the corresponding time, the proportion of damage, and the value of the affected body. Multiplying the above factors can form a dynamic risk value that is related to the probability of landslide instability, spatial exposure, temporal exposure, degree of loss, and value. In this embodiment, the dynamic risk value of a disaster-bearing body refers to the risk quantification result obtained by comprehensively considering the possibility of landslide instability, the spatial exposure degree of the disaster-bearing body, the temporal exposure degree of the disaster-bearing body, the loss rate of the disaster-bearing body, and the value of the disaster-bearing body at the time of advance simulation. For disaster-bearing bodies involving people, the dynamic risk value of the disaster-bearing body is expressed in terms of the number of people. For disaster-bearing bodies involving buildings and roads, the dynamic risk value of the disaster-bearing body is expressed in terms of monetary amount. The "dynamic" in the dynamic risk value of the disaster-bearing body means that it is updated as the rainfall intensity, water retention height, shear strength, phase shift advance time, landslide impact range, and exposure status of the disaster-bearing body change within each calculation time step. In this embodiment, at each calculation time step, the following processes are repeatedly executed: rainfall intensity input, waterlogging height update, pore water pressure update, shear strength update, strength degradation index calculation, phase shift advance time calculation, instability probability calculation, slippage influence range generation, and dynamic risk value calculation of the disaster-bearing body. , , , , , , , , , , , , and The data is synchronously stored in the digital twin contact surface model database of loess-red layer landslides, and used to update the contact surface water retention state map, strength deterioration state map, early instability probability map, sliding influence range map, and dynamic risk map of the disaster-bearing body on the display terminal. It can update the landslide instability probability, sliding influence range, and dynamic risk value of the disaster-bearing body in advance by using the strength deterioration information inside the loess-red layer contact surface before the slope surface displacement has a significant response.

[0032] The embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the present invention is not limited thereto. Various changes can be made within the scope of knowledge possessed by those skilled in the art without departing from the spirit of the present invention.

Claims

1. A method for dynamic estimation of landslide risk in loess red beds based on digital twins, characterized in that, Includes the following steps: S1. Obtain topographic profile, loess layer boundary, red layer top surface boundary, loess-red layer contact surface location, rear edge crack location, potential sliding boundary and disaster-bearing body spatial distribution data of the target landslide, construct a digital twin contact surface model of the loess-red layer landslide, and obtain basic data of the contact surface unit; S2. Input the rainfall intensity, and calculate the water retention height of the contact surface unit by combining the trailing edge crack connectivity state and the equivalent discharge coefficient of the contact surface unit in the basic data of the contact surface unit. S3. Update the pore water pressure and shear strength of the contact surface unit according to the water retention height, and calculate the strength degradation index; S4. Based on the strength degradation index and the trailing edge surface displacement, calculate the phase advance time of the strength degradation relative to the trailing edge displacement response. S5. Calculate the landslide instability probability based on the shear strength, and map the landslide instability probability to the predicted instability probability based on the phase shift lead time. S6. Generate the slippage impact range based on the instability probability and the lead time of the phase shift, and calculate the dynamic risk value of the disaster-bearing body by combining the spatial distribution data of the disaster-bearing body.

2. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 1, characterized in that, The basic data of the contact surface unit includes the contact surface unit length, contact surface inclination angle, overburden weight per unit width, initial cohesion, initial internal friction angle, initial pore water pressure, contact surface equivalent discharge coefficient, and trailing edge fracture connectivity. The trailing edge fracture connectivity includes a connected state and a disconnected state. The connected state indicates that the corresponding contact surface unit and the trailing edge fracture are hydraulically connected, while the disconnected state indicates that the corresponding contact surface unit and the trailing edge fracture are not hydraulically connected.

3. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 2, characterized in that, Calculate the water retention height of the contact surface element, including: The equivalent infiltration flux into the contact surface unit is determined based on the rainfall intensity and the connectivity of the trailing fractures. The discharge flux of the contact surface unit is determined based on the equivalent discharge coefficient of the contact surface, the water height at the previous moment, and the length of the contact surface unit. Based on the previous water level height, equivalent infiltration flux, discharge flux, and calculation time step, the current water level height is obtained according to the water conservation relationship.

4. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 3, characterized in that, Update the pore water pressure and shear strength of the contact surface elements, including: The pore water pressure at the current moment is determined based on the initial pore water pressure, the specific weight of water, and the current water level. Based on the initial cohesion, initial internal friction angle, and water content degradation relationship, determine the cohesion and internal friction angle at the current moment; The total normal stress is determined based on the component of the overburden weight per unit width in the normal direction of the contact surface and the length of the contact surface element. The shear strength at the current moment is determined according to the Mohr-Coulomb shear strength criterion based on the cohesion, total normal stress, pore water pressure, and internal friction angle at the current moment.

5. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 4, characterized in that, The strength degradation index is calculated, including: The shear strength before the start of rainfall is taken as the initial shear strength, and the ratio of the difference between the initial shear strength and the current shear strength to the initial shear strength is determined as the strength deterioration index.

6. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 5, characterized in that, Calculating the phase misalignment advance time includes: For contact surface elements in a connected state, the strength degradation index is weighted and averaged according to the length of the contact surface element to obtain the strength degradation amount of the trailing edge control area. The trailing edge displacement velocity is determined based on the difference in trailing edge surface displacement at adjacent calculation times and the calculation time step. Time-delay correlation calculations were performed on the intensity degradation of the trailing edge control zone and the trailing edge displacement velocity. The candidate advance time corresponding to the maximum correlation was determined as the phase misalignment advance time.

7. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 6, characterized in that, Calculate the probability of instability in advance, including: Based on the current shear strength and the length of the contact surface element, determine the resultant force of the anti-slip force along the contact surface according to the unit width section; The resultant force of the sliding force along the contact surface is determined based on the weight of the overlying soil per unit width and the inclination angle of the contact surface. The stability coefficient is determined based on the ratio of the resultant anti-slip force to the resultant sliding force. The probability of landslides under rainfall conditions is determined based on the correlation between the stability coefficient and the probability of landslides. The probability of occurrence of rainfall-induced factors is determined based on the frequency of occurrence of the current rainfall process or the set working conditions, and the probability of landslide instability is determined based on the probability of occurrence of rainfall-induced factors and the probability of landslide under rainfall conditions. By mapping the landslide instability probability to the projected time after the phase shift, the projected instability probability is obtained.

8. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 7, characterized in that, The generated slip influence range includes: The instability probability and the advance time of phase misalignment are input into the digital twin slip range model verified by historical numerical simulation results. The digital twin slip range model uses the maximum slip distance and the landslide instability impact range under rainfall conditions as verification objects to establish the correspondence between the instability probability and the slip impact range, and generates the slip impact range corresponding to the instability time after the advance time of phase misalignment.

9. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 8, characterized in that, Calculating the dynamic risk value of a disaster-bearing body includes: By overlaying the impact range of the slippage with the spatial distribution data of the disaster-bearing bodies, the spatial probability of disaster-bearing bodies can be determined. Determine the disaster-bearing body's temporal probability of disaster, disaster-bearing body's loss rate, and disaster-bearing body's value based on the disaster-bearing body's category; The dynamic risk value of a disaster-bearing body is determined by multiplying the instability probability, spatial disaster-bearing probability, temporal disaster-bearing probability, loss rate, and value of the disaster-bearing body in advance.

10. The method for dynamic estimation of loess red bed landslide risk based on digital twins according to claim 9, characterized in that, The disaster-bearing body categories include people, buildings, and roads; when the disaster-bearing body category is people, the disaster-bearing body value is expressed in terms of the number of people, and the dynamic risk value of the disaster-bearing body is expressed in terms of the number of people; when the disaster-bearing body category is buildings or roads, the disaster-bearing body value is expressed in terms of economic value, and the dynamic risk value of the disaster-bearing body is expressed in terms of monetary amount.

Images