Tunnel portal landslide prediction method based on shallow earthquake-stress combined inversion

Abstract

The invention belongs to the technical field of landslide prediction, and relates to a shallow earthquake-stress joint inversion tunnel portal landslide prediction method, which comprises the following steps of: obtaining slope body first-arrival wave travel time data through a high-frequency seismic source array, constructing a geologic model containing a rock mass integrity coefficient matrix, identifying rock mass structure deterioration through elastic wave velocity distribution, and predicting the landslide at a tunnel portal through shallow earthquake-stress joint inversion. The method comprises the following steps: locking a potential unstable structure weak area in advance, then embedding a rock mass integrity coefficient matrix as priori knowledge into a stress field for solving, and by establishing a coupling relationship between seismic wave propagation and stress, inverting slope stress field distribution, capturing an early stress redistribution signal caused by structural degradation of a deep rock mass, and finally obtaining an initial stress redistribution signal of the deep rock mass. Early perception of mechanical instability precursor is realized; and then calculating strain energy density distribution of the potential slip surface based on the stress and strain data, and accurately judging the potential slip surface before the slip surface is not completely cut through by judging whether the strain energy density exceeds a threshold value, whether a continuous sheet is formed and whether the expansion rate exceeds the standard, thereby avoiding lagging caused by dependence on a displacement signal.

Description

Technical Field

[0001] This application belongs to the field of landslide prediction technology, and more specifically, relates to a method for predicting landslides at tunnel entrances based on shallow earthquake-stress joint inversion. Background Technology

[0002] In mountain tunnel engineering, the delayed landslide early warning at the tunnel entrance is a core pain point restricting construction safety. This is because traditional monitoring technologies rely too heavily on displacement detection, such as GNSS surface displacement monitoring and borehole inclinometer deep displacement monitoring. These technologies use the macroscopic movement of the rock mass as the basis for judgment. However, displacement is the final manifestation of the early instability process after the deep rock mass has completed stress redistribution, structural surface failure, and slip surface connection. In order to further compress the time, higher displacement thresholds or extended observation periods are usually set. However, since these methods cannot capture the early mechanical instability signals of the gradual failure of the deep rock mass, they also lead to the problem of delayed early warning. Summary of the Invention

[0003] This invention provides a method for predicting landslides at tunnel entrances based on shallow seismic-stress joint inversion, aiming to solve the technical problem of delayed early warning caused by the inability to capture early mechanical instability signals of progressive failure of deep rock masses.

[0004] A method for predicting landslides at tunnel entrances based on shallow earthquake-stress joint inversion includes the following steps: Step S1. Deploy multiple micro-intelligent seismic source nodes in the target area to form a high-frequency seismic source array, and obtain the first arrival wave travel time data of the slope based on the high-frequency seismic source array. Construct a geological model based on the first arrival wave travel time data of the slope, and then calculate the rock mass integrity-related parameters according to the elastic wave velocity distribution to form a rock mass integrity coefficient matrix. Step S2. Fiber optic stress sensing chains are implanted in the construction borehole to capture real-time changes in the direction of principal stress and the magnitude of principal stress within the rock mass. The rock mass integrity coefficient matrix is ​​used as prior knowledge and embedded into the stress field solution process. At the same time, the correlation between seismic wave propagation and stress is established, the intrinsic relationship between stress and wave velocity is constructed, and then coupled solutions are performed to obtain the slope stress field distribution data. Step S3. Based on the stress data and strain data, calculate the strain energy density distribution at each point on the potential slip surface. Based on the calculated strain energy density, determine whether it exceeds a predetermined threshold and whether it is continuous. If it is continuous and the expansion rate of the continuous area exceeds the predetermined threshold, it is determined to be a potential slip surface. Step S4. Obtain the strain energy density exceedance rate, shear stress change rate over time, and daily propagation rate of the sliding surface. Based on the obtained indicators, perform a weighted summation to obtain the instability dynamic index, which is used for risk assessment.

[0005] This invention acquires first-arrival wave travel time data of a slope using a high-frequency seismic source array, constructs a geological model including a rock mass integrity coefficient matrix, and can identify rock mass structural deterioration through elastic wave velocity distribution before deep rock mass displacement is observed, thus identifying potential unstable structural weak areas in advance and laying a structural foundation for early signal capture. Then, the rock mass integrity coefficient matrix is ​​embedded as prior knowledge into the stress field solution, overcoming the limitations of traditional stress monitoring and isolated data interpretation. By establishing the coupling relationship between seismic wave propagation and stress, the stress field distribution of the slope is inverted, capturing early stress redistribution signals caused by structural deterioration in deep rock mass, achieving early perception of precursors to mechanical instability. Furthermore, based on stress and strain data, the strain energy density distribution of potential slip surfaces is calculated. By judging whether the strain energy density exceeds a threshold, whether it is continuous, and whether the expansion rate exceeds the standard, potential slip surfaces can be accurately identified before they are fully connected and without obvious displacement, avoiding lag caused by relying on displacement signals.

[0006] Preferably, environmental sensors and triaxial accelerometers are arranged at the multiple miniature intelligent seismic source nodes to collect rainfall intensity and construction vibration data; Adjusting the source excitation frequency based on collected rainfall intensity and construction vibration data includes: If the real-time rainfall intensity is greater than the predetermined threshold, the source excitation frequency is the product of the rainfall-frequency adjustment coefficient and the rainfall intensity plus the fundamental frequency. If the construction vibration intensity is greater than the predetermined threshold, the excitation frequency of the vibration source is the product of the vibration-frequency adjustment coefficient and the vibration intensity plus the foundation frequency. If neither of the above two situations applies, then the excitation frequency of the seismic source is the set fundamental frequency.

[0007] Preferably, the formation of the rock mass integrity coefficient matrix includes the following steps: The first arrival travel time of the slope is extracted from the received signal using an AIC automatic pickup device. Based on the extracted travel time, lithological parameters, preset weight matrix and regularization coefficient, an objective function consisting of data fitting term and model smoothing term is constructed. Then, the slowness model is corrected by combining the fast travel method with anisotropic parameters, and the theoretical forward travel time is calculated. The three-dimensional grid P-wave velocity data is output through iterative solution. When the residual exceeds the limit during the iterative solution process, the LSQR algorithm is used to update the model and introduce anisotropic constraints until the residual converges. Based on the obtained three-dimensional mesh longitudinal wave data, experimental wave velocity of intact rock mass, and rock mass structure classification criteria, the square of the ratio of inverted wave velocity to intact wave velocity is calculated for each mesh point to generate a rock mass integrity coefficient matrix; then, the rock mass state of the mesh is classified according to the numerical range of the integrity coefficient matrix, and a three-dimensional address model with rock mass state classification is constructed by combining the three-dimensional mesh.

[0008] Preferably, after the AIC automatic pickup extracts the first arrival wave travel time of the slope, it needs to correct the extracted first arrival wave travel time, including the following steps: Calculate the difference between real-time temperature and standard temperature, and the difference between real-time humidity and standard humidity. Multiply the difference in real-time temperature by the temperature calibration coefficient, add the result of the difference in humidity by the humidity calibration coefficient, and add 1 to obtain the correction factor. Multiply the correction factor by the first arrival wave travel time of the slope to obtain the corrected first arrival wave travel time of the slope.

[0009] Preferably, the intrinsic relationship between wave speed and wave speed is represented by the constructed wave speed-stress coupling equation set, which includes motion equilibrium equations, constitutive equations constrained by integrity coefficient matrix, and geometric equations. The kinematic equilibrium equations are used to describe the equilibrium relationship between rock mass density, second-order time derivative of displacement, stress divergence, and body forces. The constitutive equation constrained by the integrity coefficient matrix is ​​used to correlate stress and strain, and the integrity coefficient matrix is ​​introduced to make the elastic stiffness tensor change dynamically with the integrity coefficient matrix. The geometric equations define strain and establish the geometric relationship between displacement and strain through displacement gradient.

[0010] Preferably, the elastic stiffness tensor is corrected using an integrity coefficient: The elastic stiffness tensor is obtained by multiplying the integrity coefficient by the rock mass stiffness, adding the value of the damaged rock mass stiffness multiplied by 1 and then subtracting the integrity coefficient.

[0011] Preferably, the coupled solution to obtain the slope stress field distribution data includes the following steps: By constructing an objective function that includes stress matching terms and monitoring data constraints, and based on the constructed wave velocity-stress coupling equations, a regularized least squares optimization algorithm is used to solve the objective function until the objective function is minimized, thus obtaining a stable stress field distribution. The stress matching term is the sum of squares of the difference between the stress output by the finite element model and the measured value of the borehole gravity; the detection data constraint term is the sum of squares of the difference between the inverted stress and the FBG monitored stress, multiplied by a regularization weight.

[0012] Preferably, step S3 includes the following steps: Based on the structural plane basic parameters and the rock mass integrity coefficient matrix, multiple candidate slip surfaces are divided. The normal shear stress ratio is calculated for each candidate slip surface. If the shear stress ratio is greater than or equal to the set threshold, the corresponding candidate surface is taken as the preliminary potential slip surface. Multiply the stress tensor and strain tensor components pairwise, sum all the product results, and finally multiply the sum by a predetermined coefficient to obtain the strain energy density of the corresponding spatial point at the current time. Based on this, traverse all spatial points to generate three-dimensional strain energy density distribution data. The energy components along the gradient direction of the slip surface normal vector in the strain energy density are filtered by the Dirac function. Then, the filtered energy components are integrated over the micro-area of ​​the initial potential slip surface. The energy values ​​of all micro-areas are summed to obtain the 2D shear strain energy density distribution of the initial potential slip surface. The area of ​​the region on each slip surface with a shear strain energy density greater than or equal to the dynamic critical energy threshold is counted. The ratio of the region area to the total area of ​​the slip surface is calculated. If the ratio is greater than or equal to a predetermined ratio threshold, the energy accumulation condition is met. Then, the rate of change of the region area ratio over time is calculated. If the rate of change is greater than a predetermined rate threshold, the accelerated expansion condition is met. When both the energy accumulation condition and the accelerated expansion condition are met, the corresponding preliminary potential slip surface is determined to have entered the critical instability state and is output as the determined potential slip surface. If either of the above conditions is not met, it is determined that the critical state has not yet been entered.

[0013] Preferably, the dynamic critical energy threshold is adjusted based on the following steps: The dynamic critical energy threshold is adjusted using an adjustment function that combines structural damage correction and rainfall correction terms. The structural damage correction term is obtained by taking the natural logarithm of the ratio of the time-varying decrease of the rock mass integrity coefficient to the initial rock mass integrity coefficient, adding 1, and then multiplying it by the damage sensitivity coefficient to obtain the threshold correction ratio caused by rock mass damage. The rainfall correction term adjusts the threshold based on the cumulative rainfall using a hyperbolic tangent function; The dynamic critical energy threshold for each potential slip surface at the current time is obtained by multiplying the base threshold by 1, adding the structural damage correction term, and then multiplying by the rainfall correction term.

[0014] Preferably, the strain energy density exceedance rate is based on the actual strain energy density of the potential sliding surface minus the strain energy critical threshold to obtain the strain energy difference, and then the difference is divided by the strain energy critical threshold to obtain the strain energy density exceedance rate. The strain energy critical threshold is obtained based on the following steps: The deviation between the current rock mass integrity and the regional average level is obtained by subtracting the regional rock mass integrity benchmark value from the real-time rock mass integrity coefficient. Multiply the rock mass integrity deviation value by the correction factor and add 1 to obtain the threshold correction factor. Multiply the benchmark strain energy threshold by the threshold correction factor to obtain the corrected strain energy critical threshold.

[0015] The beneficial effects of this invention include: This invention acquires first-arrival wave travel time data of a slope using a high-frequency seismic source array, constructs a geological model including a rock mass integrity coefficient matrix, and can identify rock mass structural deterioration through elastic wave velocity distribution before deep rock mass displacement is observed, thus identifying potential unstable structural weak areas in advance and laying a structural foundation for early signal capture. Then, the rock mass integrity coefficient matrix is ​​embedded as prior knowledge into the stress field solution, overcoming the limitations of traditional stress monitoring and isolated data interpretation. By establishing the coupling relationship between seismic wave propagation and stress, the stress field distribution of the slope is inverted, capturing early stress redistribution signals caused by structural deterioration in deep rock mass, achieving early perception of precursors to mechanical instability. Furthermore, based on stress and strain data, the strain energy density distribution of potential slip surfaces is calculated. By judging whether the strain energy density exceeds a threshold, whether it is continuous, and whether the expansion rate exceeds the standard, potential slip surfaces can be accurately identified before they are fully connected and without obvious displacement, avoiding lag caused by relying on displacement signals. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of this application, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 This is an overall step diagram provided for an embodiment of the present invention.

[0018] Figure 2 The flowchart of step S1 provided in the embodiment of the present invention is shown.

[0019] Figure 3 The flowchart of step S2 provided in the embodiment of the present invention is shown.

[0020] Figure 4 The flowchart for step S3 provided in the embodiment of the present invention is shown. Detailed Implementation

[0021] To make the technical problems, technical solutions, and beneficial effects to be solved by this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and are not intended to limit the scope of this application.

[0022] See Figure 1 As shown, the method for predicting landslides at tunnel entrances based on shallow earthquake-stress joint inversion includes the following steps: Step S1. Deploy multiple micro-intelligent seismic source nodes in the target area to form a high-frequency seismic source array, and obtain the first arrival wave travel time data of the slope based on the high-frequency seismic source array. Construct a geological model based on the first arrival wave travel time data of the slope, and then calculate the rock mass integrity-related parameters according to the elastic wave velocity distribution to form a rock mass integrity coefficient matrix. See Figure 2 As shown, in this embodiment, miniature intelligent seismic source nodes are arranged at intervals of no more than 20 meters on the tunnel entrance slope and in the potential slip zone. The miniature intelligent seismic source nodes are arranged in a grid to form a three-dimensional array covering the entire detection area. Each seismic source node integrates a vibration detection sensor (such as a piezoelectric ceramic seismic source) and an environmental sensor. The environmental sensors include a sensor for acquiring rainfall intensity and a nodal triaxial accelerometer for acquiring construction vibration intensity; The working mode of the model is determined based on the obtained rainfall intensity and construction vibration intensity. The excitation frequency of the seismic source is dynamically adjusted based on the determined working mode. The specific expression is as follows: ; In the formula: Indicates the fundamental frequency; The rainfall intensity threshold is determined based on engineering experience; Indicates real-time rainfall intensity; This represents the rainfall-frequency adjustment coefficient; This represents the vibration-frequency adjustment coefficient; Indicates the intensity of construction vibration; Indicates the vibration acceleration threshold; Based on the above three working modes, we can distinguish between the rainy season mode, the construction disturbance mode, and the dry season mode, that is, when The time indicates that the current period is the rainy season; when This indicates that the current mode is construction disturbance mode; if neither exceeds the threshold, the mode will switch to dry season mode.

[0023] Furthermore, the original pulse signal of the seismic source node is modulated using a 7-bit Barker code to generate coded pulse information with anti-interference characteristics. Then, after the receiver receives the original signal containing environmental noise, it performs cross-correlation operation between the original signal and the delayed Barker code reference signal to filter out environmental noise and extract the effective signal that matches the characteristics of the Barker code reference signal. Therefore, in this embodiment, the integration of dynamic frequency modulation and coded pulse technology solves the problem of insufficient signal penetration and anti-interference ability of the seismic source in complex environments, and provides basic data support for high-quality rock mass integrity coefficient matrix for subsequent joint inversion.

[0024] The first arrival time of the slope is extracted from the received signal using an AIC automatic pickup unit. This involves dividing the received signal sequence into two segments based on the position difference of different dividing points, calculating the variance of each segment, calculating the Akaike information criterion value for each dividing point, and finally iterating through all dividing points to find the segment with the smallest Akaike information criterion value. This value represents the arrival time of the first arrival at the receiver, thus enabling automatic identification of the first arrival. The expression for calculating the Akaike information criterion value is as follows: ; In the formula: This represents the Akaike information criterion value at signal segmentation point k; Indicates the location of the signal split point; This represents the received signal sequence; Indicates the total length of the signal; Indicates the variance of the signal segment; Since temperature and humidity can interfere with the first arrival travel time data, the extracted first arrival travel time of the slope needs to be corrected based on temperature and humidity data after extraction using AIC. The specific expression is as follows: ; In the formula: This indicates the initial arrival wave travel time of the slope before correction; This indicates the travel time of the first arrival wave on the corrected slope. Indicates the temperature correction factor; Indicates the humidity correction factor; Indicates the current temperature value; Indicates the calibration ambient reference temperature; Indicates the current humidity value; This indicates the reference humidity of the calibration environment.

[0025] Based on the corrected travel time, lithological parameters, preset weight matrix, and regularization coefficients, an objective function consisting of a data fitting term and a model smoothing term is constructed: ; In the formula: This represents the slowness model vector, which is the reciprocal of the P-wave velocity. , This represents the slowness value of the i-th subsurface grid. This represents the longitudinal wave velocity of the i-th grid. This represents the total number of grids used to divide the underground space; Represents the data weight matrix; Represents the forward travel time vector; This indicates the travel time of the first arrival wave on the corrected slope. Represents the regularization coefficient; Represents the Laplace smoothing operator; Represent the objective function; Solve the equation using the fast method: ; In the formula: Represents the gradient operator; Represents the travel time scalar field; This represents the anisotropic slowness function, with spatial points as input. and direction parameters and The output is the slowness of the corresponding spatial point along a specific direction; Indicates the angle between the ray path and the normal to the fracture surface; Indicates the azimuth angle of the ray path within the fracture surface; Indicates the reference slowness; This represents the Thomsen anisotropy parameter, used to quantify the anisotropy intensity of the inclined defense line. When it is greater than 0, it indicates that the slowness increases more significantly when the ray is at a certain angle (not 0° or 90°) to the fracture normal.

[0026] Based on the objective function and equation determined above, an iterative solution is performed, with the specific steps as follows: First, an initial slow model is constructed using the assumption of homogeneity and isotropy, denoted as . The length of the vector is the same as the total number of underground grid cells. The baseline slowness is the reciprocal of the wave velocity in intact rock mass, determined based on the rock mass type in the engineering area. For example, the longitudinal wave velocity of granite in intact rock mass is fixed at 4500 m / s in the laboratory. To ensure that the initial model conforms to the physical properties of the intact rock mass; The underground three-dimensional space to be explored is uniformly divided into voxel sizes of 0.5m³ to form a discrete spatial grid system. The travel time of the spatial point where the source is located is set to 0, which means that the wave travels for 0 hours when it starts from the source. Based on this, wavefront propagation iteration is performed: the grids directly adjacent to the source point are added to the narrow band, where the narrow band is a set of grids that are temporarily stored for the travel time to be calculated; Then perform the following iterative steps until the narrowband is empty: Select the grid point with the smallest travel time in the narrow band, which is P, to represent the current position of the leading edge of the wavefront propagation; remove point P from the narrow band and add it to the receiver set, where the grid travel time in the receiver set is the determined optimal value and will not participate in subsequent updates; Calculate the travel time of all neighboring points Q of update point P: ; In the formula: This indicates the grid step size, which is 0.5m. Denotes the anisotropic slowness at the four Q points, where It can be directly calculated from the spatial directions of points P and Q; Indicate the travel time of point P; Indicates the travel time of point Q; When the narrowband is empty, the travel time calculation for all grids is completed. The travel time data of each receiver's grid is extracted to form the theoretical travel time vector. Theoretical travel vector With the corrected One-to-one correspondence; Slow-speed model correction: Step 1: Using the slowness model of the kth round as input, calculate the corresponding theoretical travel time based on the wavefront propagation iterative steps (Fast Progress Method, FMM) mentioned above. ; Step 2: Calculate the residual vector between the theoretical travel time and the measured travel time to quantify the deviation between the current model and the actual data; Step 3: If the absolute value of the residual vector is greater than the preset threshold, proceed to the next step of gradient calculation and model update. If the absolute value of the residual vector is less than or equal to a preset threshold, it indicates that the model has converged, and the slow model will be switched to a slower model. Convert to wave velocity distribution and output the wave velocity results; the threshold set here is 0.1ms. Step 4: Calculate the gradient of the objective function using the adjoint state method, where the gradient expression is as follows: ; In the formula: This represents the gradient of the objective function; The Jacobian matrix representing travel time versus slowness; This represents the transpose of the data weight matrix; Transpose of the Laplace smoothing operator; Step 5: Use the LSQR algorithm (an efficient algorithm for solving linear systems with high coefficients) to solve for the slow-speed model update. : ; Based on this updated [number] Wheel model is , This represents the slowness model for the k-th round; After each model update, the Thomsen anisotropy parameters are dynamically adjusted. Specifically, adjustments are made by analyzing the residual differences in travel time at different azimuth angles to determine the direction of parameter adjustment, while simultaneously applying physical constraints. That is, the physical value range that conforms to the anisotropy of the rock mass is returned to the first step, and the theoretical travel time is recalculated with the updated model. The next round of iteration is carried out. After 6 rounds of iteration, the residual meets the threshold requirement and the slow speed model converges. Since the slowness model is the reciprocal of the wave velocity of intact rock mass, the wave velocity distribution is obtained by iterative modification of the slowness model. The obtained wave velocity distribution is then converted into rock mass integrity coefficients, as follows: For each 3D grid, the rock mass integrity coefficient is obtained using the following formula: ; In the formula: For the inversion wave velocity corresponding to the three-dimensional grid, where Indices representing the three-dimensional mesh; Indicates the reference wave velocity of the rock mass; based on The rock mass is divided into four categories; A value greater than or equal to 0.9 indicates a complete rock mass; Rock masses with a strength greater than or equal to 0.75 and less than 0.9 are classified as slightly jointed rock masses. A rock mass with a strength greater than or equal to 0.6 and less than 0.75 is classified as a rock mass with moderately developed fractures. A value less than 0.6 is considered a fracture zone; The divided geological units are mapped to a three-dimensional spatial grid to generate a geological model in VTK format.

[0027] In this embodiment, by generating an integrity coefficient corresponding to a three-dimensional mesh, the spatial distribution of fracture zones and medium-sized fracture zones can be visually and explicitly displayed, enabling the identification of unstable and weak areas. Furthermore, rock mass deterioration can be identified through the distribution of elastic wave velocity without waiting for macroscopic displacement, thus preparing for the early identification of weak areas.

[0028] Step S2. See Figure 3 As shown, a distributed fiber optic sensor chain (FBG array) is selected as the monitoring device. Boreholes are laid out along the potential slip direction, with a borehole depth not less than the estimated depth of the sliding bed. Fiber grating sensor chains are implanted in the boreholes at a vertical spacing of 2m and a horizontal spacing of 50m. The sensor chains are spatially registered with the three-dimensional wave velocity grid in step S1. The sensor chains collect parameters such as the principal stress direction, principal stress level, and shear stress of the rock mass in real time, and integrate the collected dispersed stress data into spatially continuous stress tensor field data.

[0029] Based on the rock mass integrity coefficient and the principles of elastic dynamics and damage mechanics, three sets of interrelated coupled control equations are constructed: ; In the formula: Indicates the density of the rock mass. , Indicates the actual wave velocity of the rock mass; Represents the displacement vector field; Represents the stress tensor; Represents volume force; Represents the fourth-order elastic stiffness tensor, based on Dynamic correction; Represents the strain tensor; Represents a time variable; Represents the gradient operator; Represents the divergence operator; Indicates transpose; based on The above is corrected by using a weighted fusion method. : ; In the formula: Represents the stiffness tensor of the complete rock mass, calibrated based on indoor rock sample mechanics experiments; This represents the residual stiffness tensor of the damaged rock mass, determined based on experiments with fractured rock. Indicates the rock mass integrity coefficient; In this embodiment, The larger the value, the higher the proportion of the stiffness of the intact rock mass in the final tensor; The smaller the value, the higher the proportion of stiffness in the damaged rock mass, ensuring that the stiffness tensor can truly reflect the mechanical properties of rock masses with different integritys, enabling accurate capture of early stress concentration in deep rock masses; and through the coupling relationship between wave velocity and stress, early signals of structural deterioration to gravity distribution can be identified without waiting for displacement to manifest, solving the lag problem of traditional displacement-dependent methods.

[0030] Then, construct the objective function, and solve for the displacement vector field by minimizing the objective function. The objective function is as follows: ; In the formula: This represents the measured value of borehole stress; Represents the stiffness matrix; Indicates the regularization weight; This represents the stress tensor directly obtained from FBG monitoring; For the stress tensor that needs to be inverted; In this embodiment, the displacement vector field is used. As intermediate variables, the solution is obtained through regularized least squares optimization. Then, strain is derived through geometric equations and constitutive equations. and stress .

[0031] It should be noted that the massive amount of data from high-frequency seismic sources and FBG can easily lead to overfitting. Therefore, regularized least squares method is used to solve the problem, which balances data fitting and model stability.

[0032] Step S3. Participate Figure 4As shown, multiple candidate slip surfaces (such as areas with dense joints and obvious weathering) are divided based on the basic parameters of the structural surface and the rock mass integrity coefficient matrix. The normal shear stress ratio is calculated for each candidate slip surface. If the shear stress ratio is greater than or equal to the set threshold (0.3), the corresponding candidate surface is taken as the preliminary potential slip surface. Three-dimensional strain energy density calculation based on strain tensor field and stress tensor field: ; In the formula: This represents the stress tensor, which is the result obtained from the inversion in step S2. Here is the strain tensor, and here is the result obtained from the inversion in step S2. Represents spacetime coordinates, Indicates the location of underground space. Indicates time; For each initial potential slip surface, extract the 2D shear strain energy density distribution: ; In the formula: Represents the k-th preliminary potential slip surface 2D shear strain energy density; Represents the three-dimensional strain energy density; This represents the Dirac function; Indicates slip surface The unit normal vector; The gradient vector representing the strain energy density; Represents the area of ​​a small element of the slip surface; The area of ​​the region on each slip surface with a shear strain energy density greater than or equal to the dynamic critical energy threshold is counted. The ratio of the area of ​​the region to the total area of ​​the slip surface is calculated. If the ratio is greater than or equal to a predetermined ratio threshold (60%), the energy accumulation condition is met. The rate of change of the ratio of the area over time is then calculated. If the rate of change is greater than a predetermined rate of change threshold (0.05 / day), the accelerated expansion condition is met. When both the energy accumulation condition and the accelerated expansion condition are met, the corresponding preliminary potential slip surface is determined to have entered the critical instability state and is output as the determined potential slip surface. If either of the above conditions is not met, it is determined that the critical state has not yet been entered.

[0033] The adjustment formula for the dynamic critical energy threshold is as follows: ; ; ; In the formula: Represents the k-th slip surface The dynamic critical energy threshold at time t; Indicates the basic critical threshold; Indicates the damage sensitivity coefficient; Natural logarithm operator; Indicates slip surface nearby The time-varying decrease; This represents the initial integrity coefficient of the rock mass (benchmark value). Represents the rainfall correction function; Indicates cumulative rainfall; Indicates the initial time. Rock mass integrity coefficient at that time; This represents the rock mass integrity coefficient at time t.

[0034] In this embodiment, the structural damage correction term increases with the increase of damage, thereby reducing the dynamic critical energy threshold; the rainfall correction term increases with the increase of cumulative rainfall, thereby reducing the dynamic critical energy threshold. At the same time, the area ratio and expansion rate are added to finally identify high-risk areas that are continuously and rapidly expanding, thus avoiding interference from isolated points.

[0035] Step S4. Obtain the excess rate of strain energy density on the sliding surface. Rate of change of shear stress over time and daily expansion rate indicators

[0036] The instability dynamics index is obtained by weighted summation based on the acquired indicators. The specific expression is as follows: ; In the formula: This represents the difference between the actual strain energy density of the sliding surface and the dynamic critical energy threshold. This represents the maximum shear stress on the potential slip surface (real-time monitoring by fiber optic stress sensing chain). Represents the continuous area of ​​the region where strain energy exceeds the limit; , as well as Indicates the weighting coefficient; The dynamic critical energy threshold; based on Different early warning mechanisms can be set, for example A yellow alert will be issued if the value is between 0.3 and 0.6. If the value is in the range of 0.6-0.8, an orange alert will be issued; if... A red alert is issued if the value is greater than 0.8. Different response measures are set for different alert levels. The setting of response measures is a conventional technical means in this field, so it will not be described in detail in this invention.

[0037] The above are merely preferred embodiments of this application and are not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for predicting landslides at tunnel entrances based on shallow seismic-stress joint inversion, characterized in that, Includes the following steps: Step S1. Deploy multiple micro-intelligent seismic source nodes in the target area to form a high-frequency seismic source array, and obtain the first arrival wave travel time data of the slope based on the high-frequency seismic source array. Construct a geological model based on the first arrival wave travel time data of the slope, and then calculate the rock mass integrity-related parameters according to the elastic wave velocity distribution to form a rock mass integrity coefficient matrix. Step S2. Fiber optic stress sensing chains are implanted in the construction borehole to capture real-time changes in the direction of principal stress and the magnitude of principal stress within the rock mass. The rock mass integrity coefficient matrix is ​​used as prior knowledge and embedded into the stress field solution process. At the same time, the correlation between seismic wave propagation and stress is established, the intrinsic relationship between stress and wave velocity is constructed, and then coupled solutions are performed to obtain the slope stress field distribution data. Step S3. Based on the stress data and strain data, calculate the strain energy density distribution at each point on the potential slip surface. Based on the calculated strain energy density, determine whether it exceeds a predetermined threshold and whether it is continuous. If it is continuous and the expansion rate of the continuous area exceeds the predetermined threshold, it is determined to be a potential slip surface. Step S4. Obtain the strain energy density exceedance rate, shear stress change rate over time, and daily propagation rate of the sliding surface. Based on the obtained indicators, perform a weighted summation to obtain the instability dynamic index, which is used for risk assessment.

2. The method for predicting tunnel entrance landslides based on shallow seismic-stress joint inversion according to claim 1, characterized in that, Environmental sensors and triaxial accelerometers were deployed at the multiple miniature intelligent seismic source nodes to collect rainfall intensity and construction vibration data. Adjusting the source excitation frequency based on collected rainfall intensity and construction vibration data includes: If the real-time rainfall intensity is greater than the predetermined threshold, the source excitation frequency is the product of the rainfall-frequency adjustment coefficient and the rainfall intensity plus the fundamental frequency. If the construction vibration intensity is greater than the predetermined threshold, the excitation frequency of the vibration source is the product of the vibration-frequency adjustment coefficient and the vibration intensity plus the foundation frequency. If neither of the above two situations applies, then the excitation frequency of the seismic source is the set fundamental frequency.

3. The method for predicting tunnel entrance landslides using shallow seismic-stress joint inversion according to claim 1, characterized in that, The formation of the rock mass integrity coefficient matrix includes the following steps: The first arrival travel time of the slope is extracted from the received signal using an AIC automatic pickup device. Based on the extracted travel time, lithological parameters, preset weight matrix and regularization coefficient, an objective function consisting of data fitting term and model smoothing term is constructed. Then, the slowness model is corrected by combining the fast travel method with anisotropic parameters, and the theoretical forward travel time is calculated. The three-dimensional grid P-wave velocity data is output through iterative solution. When the residual exceeds the limit during the iterative solution process, the LSQR algorithm is used to update the model and introduce anisotropic constraints until the residual converges. Based on the obtained three-dimensional mesh longitudinal wave data, experimental wave velocity of intact rock mass, and rock mass structure classification criteria, the square of the ratio of inverted wave velocity to intact wave velocity is calculated for each mesh point to generate a rock mass integrity coefficient matrix; then, the rock mass state of the mesh is classified according to the numerical range of the integrity coefficient matrix, and a three-dimensional address model with rock mass state classification is constructed by combining the three-dimensional mesh.

4. The method for predicting tunnel entrance landslides using shallow seismic-stress joint inversion according to claim 3, characterized in that, After the AIC automatic pickup extracts the first arrival wave travel time of the slope, the extracted first arrival wave travel time needs to be corrected, including the following steps: Calculate the difference between real-time temperature and standard temperature, and the difference between real-time humidity and standard humidity. Multiply the difference in real-time temperature by the temperature calibration coefficient, add the result of the difference in humidity by the humidity calibration coefficient, and add 1 to obtain the correction factor. Multiply the correction factor by the first arrival wave travel time of the slope to obtain the corrected first arrival wave travel time of the slope.

5. The method for predicting tunnel entrance landslides using shallow seismic-stress joint inversion according to claim 1, characterized in that, The intrinsic relationship between wave speeds is expressed based on the constructed wave speed-stress coupling equation set, which includes motion equilibrium equations, constitutive equations constrained by integrity coefficient matrix, and geometric equations. The kinematic equilibrium equations are used to describe the equilibrium relationship between rock mass density, second-order time derivative of displacement, stress divergence, and body forces. The constitutive equation constrained by the integrity coefficient matrix is ​​used to correlate stress and strain, and the integrity coefficient matrix is ​​introduced to make the elastic stiffness tensor change dynamically with the integrity coefficient matrix. The geometric equations define strain and establish the geometric relationship between displacement and strain through displacement gradient.

6. The method for predicting tunnel entrance landslides based on shallow seismic-stress joint inversion according to claim 5, characterized in that, For the elastic stiffness tensor, an integrity coefficient is used for correction: The elastic stiffness tensor is obtained by multiplying the integrity coefficient by the rock mass stiffness, adding the value of the damaged rock mass stiffness multiplied by 1 and then subtracting the integrity coefficient.

7. The method for predicting tunnel entrance landslides based on shallow seismic-stress joint inversion according to claim 5, characterized in that, The coupled solution to obtain the slope stress field distribution data includes the following steps: By constructing an objective function that includes stress matching terms and monitoring data constraints, and based on the constructed wave velocity-stress coupling equations, a regularized least squares optimization algorithm is used to solve the objective function until the objective function is minimized, thus obtaining a stable stress field distribution. The stress matching term is the sum of squares of the difference between the stress output by the finite element model and the measured value of the borehole gravity; the detection data constraint term is the sum of squares of the difference between the inverted stress and the FBG monitored stress, multiplied by a regularization weight.

8. The method for predicting tunnel entrance landslides based on shallow seismic-stress joint inversion according to claim 1, characterized in that, Step S3 includes the following steps: Based on the structural plane basic parameters and the rock mass integrity coefficient matrix, multiple candidate slip surfaces are divided. The normal shear stress ratio is calculated for each candidate slip surface. If the shear stress ratio is greater than or equal to the set threshold, the corresponding candidate surface is taken as the preliminary potential slip surface. Multiply the stress tensor and strain tensor components pairwise, sum all the product results, and finally multiply the sum by a predetermined coefficient to obtain the strain energy density of the corresponding spatial point at the current time. Based on this, traverse all spatial points to generate three-dimensional strain energy density distribution data. The energy components along the gradient direction of the slip surface normal vector in the strain energy density are filtered by the Dirac function. Then, the filtered energy components are integrated over the micro-area of ​​the initial potential slip surface. The energy values ​​of all micro-areas are summed to obtain the 2D shear strain energy density distribution of the initial potential slip surface. The area of ​​the region on each slip surface with a shear strain energy density greater than or equal to the dynamic critical energy threshold is counted. The ratio of the region area to the total area of ​​the slip surface is calculated. If the ratio is greater than or equal to a predetermined ratio threshold, the energy accumulation condition is met. Then, the rate of change of the region area ratio over time is calculated. If the rate of change is greater than a predetermined rate threshold, the accelerated expansion condition is met. When both the energy accumulation condition and the accelerated expansion condition are met, the corresponding preliminary potential slip surface is determined to have entered the critical instability state and is output as the determined potential slip surface. If either of the above conditions is not met, it is determined that the critical state has not yet been entered.

9. The method for predicting tunnel entrance landslides based on shallow seismic-stress joint inversion according to claim 8, characterized in that, The dynamic critical energy threshold is adjusted based on the following steps: The dynamic critical energy threshold is adjusted using an adjustment function that combines structural damage correction and rainfall correction terms. The structural damage correction term is obtained by taking the natural logarithm of the ratio of the time-varying decrease of the rock mass integrity coefficient to the initial rock mass integrity coefficient, adding 1, and then multiplying it by the damage sensitivity coefficient to obtain the threshold correction ratio caused by rock mass damage. The rainfall correction term adjusts the threshold based on the cumulative rainfall using a hyperbolic tangent function; The dynamic critical energy threshold for each potential slip surface at the current time is obtained by multiplying the base threshold by 1, adding the structural damage correction term, and then multiplying by the rainfall correction term.

10. The method for predicting tunnel entrance landslides based on shallow seismic-stress joint inversion according to claim 1, characterized in that, The strain energy density exceedance rate is based on the actual strain energy density of the potential sliding surface minus the strain energy critical threshold to obtain the strain energy difference, and then the difference is divided by the strain energy critical threshold to obtain the strain energy density exceedance rate. The strain energy critical threshold is obtained based on the following steps: The deviation between the current rock mass integrity and the regional average level is obtained by subtracting the regional rock mass integrity benchmark value from the real-time rock mass integrity coefficient. Multiply the rock mass integrity deviation value by the correction factor and add 1 to obtain the threshold correction factor. Multiply the benchmark strain energy threshold by the threshold correction factor to obtain the corrected strain energy critical threshold.

Images