A method for predicting landslides at tunnel entrances using shallow earthquake-stress joint inversion.

Abstract

This application belongs to the field of landslide prediction technology, and relates to a method for predicting landslides at tunnel entrances using shallow seismic-stress joint inversion. This invention acquires first-arrival wave travel time data of the slope using a high-frequency seismic source array, constructs a geological model including a rock mass integrity coefficient matrix, identifies rock mass structural deterioration through elastic wave velocity distribution, and identifies potential unstable structural weak areas in advance. Then, the rock mass integrity coefficient matrix is ​​embedded as prior knowledge into the stress field solution. 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 masses, achieving early detection of precursors to mechanical instability. Furthermore, based on stress and strain data, the strain energy density distribution of the potential slip surface is calculated. By judging whether the strain energy density exceeds a threshold, whether it is continuous, and whether the expansion rate exceeds the standard, the potential slip surface can be accurately identified before it is fully connected, avoiding the lag caused by relying on displacement signals.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of landslide prediction, and more particularly to a tunnel portal landslide prediction method based on shallow seismic-stress joint inversion. BACKGROUND

[0002] In mountainous tunnel engineering, the delay of tunnel portal landslide warning is the core pain point that restricts construction safety. The traditional monitoring technology excessively relies on displacement detection, such as GNSS ground surface displacement monitoring and borehole inclinometer deep displacement monitoring, which takes the macroscopic movement of rock mass as the basis for judgment. However, displacement is the final manifestation after the early instability process of deep rock mass, such as stress redistribution, structure surface damage, and slip surface penetration. In order to further compress the time, a higher displacement threshold or a longer observation period is usually set. However, due to the inability to capture the early mechanical instability signals of the progressive damage of deep rock mass, the problem of delayed warning is also caused. SUMMARY

[0003] The present application provides a tunnel portal landslide prediction method based on shallow seismic-stress joint inversion, which aims to solve the technical problem of delayed warning caused by the inability to capture the early mechanical instability signals of the progressive damage of deep rock mass.

[0004] The tunnel portal landslide prediction method based on shallow seismic-stress joint inversion comprises the following steps:

[0005] Step S1. Arranging multiple micro-intelligent seismic source nodes in the target area to form a high-frequency seismic source array, and obtaining slope first arrival wave travel time data based on the high-frequency seismic source array. Based on the slope first arrival wave travel time data, a geological model is constructed, and then the rock mass integrity related parameters are calculated according to the elastic wave velocity distribution to form a rock mass integrity coefficient matrix.

[0006] Step S2. Implanting a fiber Bragg grating stress sensing chain in the construction drill hole to capture real-time change data of the main stress direction and the main stress level of the rock mass inside. The rock mass integrity coefficient matrix is used as prior knowledge and embedded in the stress field solving process. At the same time, the relationship between seismic wave propagation and stress is established, and the internal relationship between stress and wave velocity is constructed. Then, the coupling solution is performed to obtain the slope stress field distribution data.

[0007] Step S3. According to the stress data and strain data, the strain energy density distribution of each point on the potential slip surface is calculated. Based on the calculated strain energy density, it is determined whether it exceeds the predetermined threshold and whether it is continuous and in pieces. If it is continuous and in pieces, and the expansion rate of the piece region exceeds the predetermined threshold, it is determined as a potential slip surface.

[0008] Step S4. Obtain the strain energy density exceeding rate of the sliding surface, the change rate of the shear stress with time and the daily extension rate index, and perform weighted summation based on the obtained indexes to obtain an instability dynamics index for risk assessment.

[0009] The present application can identify rock mass structure degradation through elastic wave velocity distribution before deep rock mass displacement appears, and lock potential instability structure weak area in advance, laying a structural foundation for early signal capture; then, the rock mass integrity coefficient matrix is embedded as prior knowledge in stress field solving, breaking through the limitation of traditional stress monitoring in isolated data interpretation, and through the establishment of the coupling relationship between seismic wave propagation and stress, the slope stress field distribution is inverted, the early stress redistribution signal caused by structure degradation of deep rock mass is captured, and early perception of mechanical instability precursor is realized; then, based on stress and strain data, the strain energy density distribution of potential sliding surface is calculated, and through judging whether the strain energy density exceeds the threshold, whether it is continuous and whether the expansion rate is excessive, the potential sliding surface can be accurately judged before the sliding surface is completely penetrated and obvious displacement occurs, avoiding the lag caused by relying on displacement signal.

[0010] Preferably, the multiple micro intelligent seismic source nodes are arranged with environmental sensors and triaxial accelerometers to collect rainfall intensity and construction vibration data.

[0011] The seismic source excitation frequency is adjusted based on the collected rainfall intensity and construction vibration data, including:

[0012] If the real-time rainfall intensity is greater than a predetermined threshold, the seismic source excitation frequency is the product of the rainfall-frequency adjustment coefficient and the rainfall intensity plus the basic frequency;

[0013] If the construction vibration intensity is greater than a predetermined threshold, the seismic source excitation frequency is the product of the vibration-frequency adjustment coefficient and the vibration intensity plus the basic frequency;

[0014] If neither of the above two cases occurs, the seismic source excitation frequency is the set basic frequency.

[0015] Preferably, the formation of the rock mass integrity coefficient matrix includes the following steps:

[0016] Based on the AIC automatic picker, the slope first arrival wave travel time is extracted from the received signal, based on the extracted travel time, lithology parameters, a preset weight matrix and a regularization coefficient, a target function composed of a data fitting term and a model smoothing term is constructed, then the fast marching method is combined with anisotropy parameter correction slowness model to calculate the theoretical forward travel time, and three-dimensional grid P-wave velocity data is output through iterative solving, in the iterative solving process, when the residual exceeds the limit, the LSQR algorithm is used to update the model and introduce anisotropy constraint, until the residual converges;

[0017] Based on the obtained three-dimensional grid P-wave data, complete rock mass experimental wave velocity and rock mass structure classification standard, the square of the ratio of the inversion wave velocity to the complete wave velocity is calculated for each grid point to generate a rock mass integrity coefficient matrix; then the grid rock mass state 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 in combination with the three-dimensional grid.

[0018] Preferably, after the AIC automatic picker extracts the slope first arrival wave travel time, the extracted slope first arrival wave travel time needs to be corrected, including the following steps:

[0019] The difference between the real-time temperature and the standard temperature, and the difference between the real-time humidity and the standard humidity are calculated respectively, and the correction factor is obtained by adding 1 to the result of multiplying the difference in real-time temperature by the temperature calibration coefficient and then adding the product of the humidity difference and the humidity calibration coefficient, and the corrected slope first arrival wave travel time is obtained by multiplying the correction factor by the slope first arrival wave travel time.

[0020] Preferably, the internal relationship between the wave velocities is based on the constructed wave velocity-stress coupling equation set, wherein the wave velocity-stress coupling equation set includes a motion balance equation, a constitutive equation constrained by an integrity coefficient matrix, and a geometric equation;

[0021] The motion balance equation is used to describe the balance relationship between the rock mass density, the second-order time derivative of displacement, and the stress divergence and volume force;

[0022] The constitutive equation constrained by the integrity coefficient matrix is used to relate stress and strain, and the integrity coefficient matrix is introduced to make the elastic stiffness tensor dynamically change with the integrity coefficient matrix;

[0023] The geometric equation defines the strain through the displacement gradient and establishes the geometric relationship between the displacement and the strain.

[0024] Preferably, for the elastic stiffness tensor, the integrity coefficient is used for correction:

[0025] The elastic stiffness tensor is obtained by multiplying the integrity coefficient by the rock mass stiffness and adding the value obtained by multiplying the damaged rock mass stiffness by 1 minus the integrity coefficient.

[0026] Preferably, the coupling solution to obtain the slope stress field distribution data includes the following steps:

[0027] A target function is constructed by including a stress matching term and a monitoring data constraint term, and based on the constructed wave velocity-stress coupling equation set, a regularization least squares optimization algorithm is used to solve the target function until the target function is minimized to obtain a stable stress field distribution;

[0028] Wherein the stress matching term is the sum of squares of differences between the stress calculated by the finite element model and the measured value of the drilling gravity; and the detection data constraint term is the sum of squares of differences between the inverted stress and the FBG monitoring stress, multiplied by a regularization weight.

[0029] Preferably, the step S3 comprises the following steps:

[0030] Based on the structural plane basic parameters and the rock mass integrity coefficient matrix, a plurality of candidate slip surfaces are divided, and the normal shear stress ratio is calculated for each candidate slip surface. If the shear stress ratio is greater than or equal to a set threshold value, the corresponding candidate surface is taken as a preliminary potential slip surface.

[0031] The components corresponding to the stress tensor and the strain tensor are multiplied two by two, and then all the product results are summed up. Finally, the sum result is multiplied by a predetermined coefficient to obtain the strain energy density of the corresponding spatial point at the current time. Based on this, all spatial points are traversed to generate three-dimensional strain energy density distribution data.

[0032] The energy component along the gradient direction of the normal vector of the slip surface in the strain energy density is filtered through the Dirac function, and then the filtered energy component is integrated on the microelement area of the preliminary potential slip surface. The energy values of all microelement areas are accumulated to obtain the 2D shear strain energy density distribution of the preliminary potential slip surface.

[0033] The area of the region on each slip surface where the shear strain energy density is greater than or equal to the dynamic critical energy threshold value is counted, and the ratio of the area to the total area of the slip surface is calculated. If the ratio is greater than or equal to a predetermined proportion threshold value, the energy accumulation condition is met. Then the change rate of the area ratio with time is calculated. If the change rate is greater than a predetermined rate threshold value, the accelerated expansion condition is met. When both the energy accumulation condition and the accelerated expansion condition are met, it is determined that the corresponding preliminary potential slip surface enters a critical unstable state, and it is output as a determined potential slip surface. If either of the above conditions is not met, it is determined that it has not entered a critical state.

[0034] Preferably, the dynamic critical energy threshold value is adjusted based on the following steps:

[0035] An adjustment function combining a structural damage correction term and a rainfall correction term is used to adjust the dynamic critical energy threshold value. The structural damage correction term is obtained by taking the natural logarithm of the ratio of the time-varying amplitude reduction of the rock mass integrity coefficient to the initial rock mass integrity coefficient, adding 1, and then multiplying by a damage sensitivity coefficient to obtain the threshold value correction proportion caused by rock mass damage.

[0036] The rainfall correction term adjusts the threshold value according to the cumulative rainfall amount through a hyperbolic tangent function.

[0037] The dynamic critical energy threshold of each potential sliding surface at the current time is obtained by multiplying the basic threshold value by 1 plus the sum of the structure damage correction terms, and then multiplying the rainfall correction term.

[0038] Preferably, the strain energy density overage rate is based on the actual strain energy density of the potential sliding surface minus the strain energy critical threshold value, to obtain a strain energy difference value, and then dividing the difference value by the strain energy critical threshold value to obtain the strain energy density overage rate.

[0039] The strain energy critical threshold value is obtained based on the following steps:

[0040] The deviation value of the current rock mass integrity from the regional average level is obtained based on the real-time rock mass integrity coefficient minus the regional rock mass integrity benchmark value.

[0041] The rock mass integrity deviation value is multiplied by a correction coefficient, and then 1 is added to obtain a threshold correction factor, and the benchmark strain energy threshold value is multiplied by the threshold correction factor to obtain the corrected strain energy critical threshold value.

[0042] The beneficial effects of the present application include:

[0043] The present application obtains the first arrival wave travel time data of the slope body through the high-frequency seismic source array, constructs a geological model containing the rock mass integrity coefficient matrix, can identify the rock mass structure degradation through the elastic wave velocity distribution before the deep rock mass appears displacement, and locks the potential unstable structure weak area in advance, lays the structural foundation for early signal capture; then embeds the rock mass integrity coefficient matrix into the stress field solution as prior knowledge, breaks through the limitations of traditional stress monitoring isolated data interpretation, establishes the coupling relationship between the seismic wave propagation and the stress, inverts the stress field distribution of the slope body, captures the early stress redistribution signal of the deep rock mass caused by structure degradation, and realizes the early perception of mechanical instability precursor; then calculates the strain energy density distribution of the potential sliding surface based on the stress and strain data, and through judging whether the strain energy density exceeds the threshold value, whether it is continuous and whether the expansion rate is overage, the potential sliding surface can be accurately judged before the sliding surface is completely penetrated and there is no obvious displacement, avoiding the lag caused by relying on displacement signals. BRIEF DESCRIPTION OF DRAWINGS

[0044] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0045] Figure 1 The overall step block diagram provided by the embodiments of the present application.

[0046] Figure 2The flow block diagram of step S1 provided for the embodiment of the present application is shown in the following.

[0047] Figure 3 The flow block diagram of step S2 provided for the embodiment of the present application is shown in the following.

[0048] Figure 4 The flow block diagram of step S3 provided for the embodiment of the present application is shown in the following. DETAILED DESCRIPTION

[0049] In order to make the technical problems, technical solutions and beneficial effects of the present application clearer, the present application will be further described in detail below in combination with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and not to limit the present application.

[0050] Referring to Figure 1 As shown in the figure, the tunnel portal landslide prediction method of the shallow seismic-stress joint inversion comprises the following steps:

[0051] Step S1. Arranging a plurality of micro intelligent seismic source nodes in the target area to form a high-frequency seismic source array, and acquiring slope first arrival wave travel time data based on the high-frequency seismic source array, constructing a geological model based on the slope first arrival wave travel time data, and then calculating rock mass integrity related parameters according to the elastic wave velocity distribution to form a rock mass integrity coefficient matrix;

[0052] Referring to Figure 2 As shown in the figure, in the present embodiment, the micro intelligent seismic source nodes are arranged at a spacing of not more than 20 meters on the tunnel portal slope and the potential sliding area, and the micro intelligent seismic source nodes are arranged in a grid to form a three-dimensional array covering the entire detection area; each seismic source node is internally integrated with a vibration detection sensor (such as a piezoelectric ceramic seismic source) and an environmental sensor;

[0053] The environmental sensor comprises a sensor for acquiring rainfall intensity and a node triaxial accelerometer for acquiring construction vibration intensity;

[0054] Based on the acquired rainfall intensity and construction vibration intensity, the working mode of the model is judged, and the seismic source excitation frequency is dynamically adjusted based on the judged working mode, and the specific expression is as follows: ;

[0055] In the formula: The basic frequency is represented by f0; The rainfall intensity threshold is represented by f1, which is determined based on engineering experience; The real-time rainfall intensity is represented by f2; The rainfall-frequency adjustment coefficient is represented by f3; The vibration-frequency adjustment coefficient is represented by f4; The construction vibration intensity is represented by f5; The vibration acceleration threshold is represented by f6;

[0056] Based on the above three working modes, the rainy season mode, the construction disturbance mode and the dry season mode are distinguished, that is, when , it indicates that the current is the rainy season mode; when , it indicates that the current is the construction disturbance mode; if neither of them exceeds the threshold, the dry season mode is entered.

[0057] Further, the original pulse signal of the seismic source node is modulated by using a 7-bit Barker code to generate coded pulse information with anti-interference characteristics; then, after receiving the original signal containing environmental noise, the receiving end performs cross-correlation operation on the original signal and the delayed Barker code reference signal, filters the environmental noise, and extracts the effective signal matched with the characteristics of the Barker code reference signal;

[0058] Therefore, in the embodiment, based on the fusion of dynamic frequency modulation and coded pulse technology, the deficiencies of signal penetration and anti-interference ability of the seismic source in a complex environment are solved, and high-quality rock mass integrity coefficient matrix basic data support is provided for subsequent joint inversion.

[0059] Based on the AIC automatic picker, the slope first arrival wave travel time is extracted from the received signal, that is, the received signal sequence is differentially divided into two segments according to different segmentation point positions, the variance of each segment signal is calculated, then the Akaike information criterion value is calculated for each segmentation point, and finally all segmentation points are traversed to find the segmentation point with the minimum Akaike information criterion value, which is the time when the first arrival wave reaches the receiving end. Based on this, automatic identification of the first arrival wave is realized; the calculation expression of the Akaike information criterion value is as follows:

[0060] ;

[0061] In the formula: Ak(k) represents the Akaike information criterion value of the signal segmentation point k; k represents the signal segmentation point position; s(k) represents the received signal sequence; N represents the total length of the signal; var(k) represents the signal segment variance;

[0062] Since temperature and humidity will interfere with the first arrival wave travel time data, after extracting the slope first arrival wave travel time based on AIC, the extracted first arrival wave travel time needs to be corrected based on temperature and humidity data, and the specific expression is as follows: ;

[0063] In the formula: t0 represents the slope first arrival wave travel time before correction; t represents the slope first arrival wave travel time after correction; a represents the temperature correction coefficient; b represents the humidity correction coefficient; T represents the current temperature value; Indicates the calibration ambient reference temperature; Indicates the current humidity value; This indicates the reference humidity of the calibration environment.

[0064] 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: ;

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

[0066] Solve the equation using the fast method: ;

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

[0068] Based on the objective function and equation determined above, an iterative solution is performed, with the specific steps as follows:

[0069] 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. For the reference slowness, the inverse of the wave velocity of intact rock, which is determined according to the rock type in the engineering area, such as granite is fixed at 4500m / s in the laboratory longitudinal wave velocity of intact rock, so , to ensure that the initial model meets the physical properties of intact rock;

[0070] The underground three-dimensional space to be detected is uniformly divided according to the voxel specification of 0.5m³ to form a discrete spatial grid system, and the travel time of the spatial point where the source is located is set to 0, indicating that the propagation time is 0 when the wave starts from the source;

[0071] Based on this, the wavefront propagation iteration is carried out: the grid directly adjacent to the source point is added to the narrow band, wherein the narrow band is a collection of grids to be calculated travel time;

[0072] The following iteration steps are executed again until the narrow band is empty:

[0073] Select the grid point with the smallest travel time in the narrow band, that is, P, which represents the current wavefront propagation front position; remove point P from the narrow band and add it to the acceptance set, wherein the grid travel time in the acceptance set is the optimal value determined and no longer participates in subsequent updates;

[0074] Calculate the travel time of all neighbor points Q of the update point P: ;

[0075] In the formula: represents the grid step, that is, 0.5m; represents the anisotropic slowness of the four Q points, wherein is directly calculated from the spatial direction of points P and Q; represents the travel time of point P; represents the travel time of point Q;

[0076] When the narrow band is empty, the travel time calculation of all grids is completed, and the travel time data of the grid where each receiver is located is extracted to form a theoretical travel time vector The theoretical travel time vector corresponds to the corrected one by one;

[0077] Slowness model correction:

[0078] First step: take the kth round slowness model as input, and calculate the corresponding theoretical travel time based on the above-mentioned wavefront propagation iteration step (fast marching method FMM);

[0079] Second step: calculate the residual vector of the theoretical travel time and the measured travel time, and quantify the deviation between the current model and the actual data;

[0080] Third step: if the absolute value of residual vector is greater than the preset threshold, then go to the next step of gradient calculation and model updating;

[0081] If the absolute value of residual vector is less than or equal to the preset threshold, it indicates that the model converges, and the slowness model is converted into wave velocity distribution, and the wave velocity result is output; the threshold set here is 0.1 ms;

[0082] Fourth step: the gradient of the objective function is calculated using the adjoint state method, and the gradient expression is as follows:

[0083] ;

[0084] In the formula: the gradient of the objective function is represented; the Jacobian matrix of travel time to slowness is represented; the transpose of the data weight matrix is represented; the transpose of the Laplace smoothing operator is represented;

[0085] Fifth step: the LSQR algorithm (efficient solution coefficient linear system algorithm) is used to solve the slowness model update amount :

[0086] ;

[0087] The updated model of the first round based on this is , wherein k represents the kth round of slowness model;

[0088] After each round of model updating, the Thomsen anisotropy parameter is dynamically adjusted, and the adjustment is specifically determined by analyzing the residual difference of travel time at different azimuth angles, and the parameter adjustment direction is determined, and a physical constraint is applied, that is, the physical value range of the rock mass anisotropy is met, and the model is returned to the first step to recalculate the theoretical travel time with the updated model, and the next round of iteration is performed, and after 6 rounds of iteration, the residual meets the threshold requirement, and the slowness model converges;

[0089] Since the slowness model is the inverse of the complete rock mass wave velocity, the wave velocity distribution is obtained by modifying and iterating the slowness model, and the obtained wave velocity distribution is converted into the rock mass integrity coefficient, and the specific process is as follows:

[0090] For each three-dimensional grid, the rock mass integrity coefficient is obtained using the following formula: ;

[0091] In the formula: is the inversion wave velocity corresponding to the three-dimensional grid, wherein the index of the three-dimensional grid is represented; represents the reference wave velocity of rock mass;

[0092] Based on The rock mass is divided into 4 categories; greater than or equal to 0.9 is determined as complete rock mass; greater than or equal to 0.75 and less than 0.9 is determined as slightly jointed rock mass; greater than or equal to 0.6 and less than 0.75 is determined as medium fissure developed rock mass; less than 0.6 is determined as fracture zone;

[0093] The divided geological unit is corresponded with the three-dimensional space grid to generate a geological model in VTK format.

[0094] In this embodiment, by generating the integrity coefficient corresponding to the three-dimensional grid, the spatial distribution of the fracture zone and the medium fissure zone can be directly and visually displayed, so that the unstable weak zone can be locked, and the rock mass degradation can be identified through the elastic wave velocity distribution, without waiting for the macroscopic displacement, so as to prepare for locking the weak zone in advance.

[0095] Step S2. Referring to Figure 3 As shown in the figure, a distributed optical fiber sensing chain (FBG array) is selected as the monitoring device, a borehole is arranged along the potential sliding direction, the borehole depth is not less than the estimated depth of the slide bed, and the optical fiber grating sensing chain is implanted in the borehole at a vertical interval of 2m and a horizontal interval of 50m, wherein the sensing chain is spatially registered with the three-dimensional wave velocity grid in step S1; the rock mass principal stress direction, principal stress level, shear stress and other parameters are collected in real time through the sensing chain, and the collected dispersed stress data is integrated into spatially continuous stress tensor field data.

[0096] Based on the rock mass integrity coefficient, three groups of coupled control equations are constructed based on the principles of elastic dynamics and damage mechanics:

[0097] In the formula: represents the density of the rock mass, , represents the actual wave velocity of the rock mass; represents the displacement vector field; represents the stress tensor; represents the body force; represents the fourth-order elastic stiffness tensor, which is dynamically corrected based on ; represents the strain tensor; represents the time variable; represents the gradient operator; represents the divergence operator; represents the transpose;

[0098] Based on ​The weighted fusion method is used to correct the stress tensor :

[0099] ;

[0100] In the formula: represents the complete rock mass stiffness tensor, which is calibrated based on laboratory rock sample mechanical experiments; represents the residual stiffness tensor of the damaged rock mass, which is determined based on broken rock experiments; represents the rock mass integrity coefficient;

[0101] In the embodiment, the value of the rock mass integrity coefficient is greater than 1, and the value of the rock mass integrity coefficient is less than 1. The greater the value is, the higher the proportion of the complete rock mass stiffness in the final tensor is. The smaller the value is, the higher the proportion of the damaged rock mass stiffness is, which ensures that the stiffness tensor can truly reflect the mechanical properties of different integrity rock masses, so that the early stress concentration of deep rock mass can be accurately captured; and through the coupling relationship between wave velocity and stress, early signals of structural degradation to gravity distribution are identified, without waiting for displacement to appear, solving the lag problem of traditional displacement dependence.

[0102] The target function is reconstructed, and the displacement vector field is solved by minimizing the target function , wherein the target function is as follows:

[0103] ;

[0104] In the formula: represents the measured value of the borehole stress; represents the stiffness matrix; represents the regularization weight; represents the stress tensor directly obtained by the FBG monitoring; is the stress tensor to be inverted;

[0105] In the embodiment, the displacement vector field is used as an intermediate variable, and the regularization least square optimization is used to solve , and then the strain and the stress are derived through the geometric equation and the constitutive equation.

[0106] It should be noted that due to the mass data of the high-frequency seismic source and the FBG, overfitting is easy to occur, and therefore the regularization least square method is used to solve, so as to balance data fitting and model stability.

[0107] Step S3. Participate Figure 4As shown, based on the structural plane basic parameters and the rock mass integrity coefficient matrix, multiple candidate slip surfaces (such as areas with dense joints and obvious weathering) are divided, the normal shear stress ratio is calculated for each candidate slip surface, and if the shear stress ratio is greater than or equal to a set threshold value (0.3), the corresponding candidate surface is taken as a preliminary potential slip surface;

[0108] Based on the strain tensor field and the stress tensor field, three-dimensional strain energy density calculation is performed: ;

[0109] In the formula: denotes the stress tensor, which is the result obtained in step S2; denotes the strain tensor, which is the result obtained in step S2; denotes the space-time coordinates, denotes the underground space position, denotes time;

[0110] For each preliminary potential slip surface, 2D shear strain energy density distribution is extracted: ;

[0111] In the formula: denotes the 2D shear strain energy density of the kth preliminary potential slip surface ; denotes the three-dimensional strain energy density; denotes the Dirac function; denotes the unit normal vector of the slip surface ; denotes the gradient vector of the strain energy density; denotes the infinitesimal area of the slip surface;

[0112] The area of the region on each slip surface where the shear strain energy density is greater than or equal to the dynamic critical energy threshold value is counted, the ratio of the area to the total area of the slip surface is calculated, and if the ratio is greater than or equal to a predetermined proportion threshold value (60%), the energy accumulation condition is met; the rate of change of the area ratio with time is calculated, and if the rate of change is greater than a predetermined change rate threshold value (0.05 / day), the accelerated expansion condition is met; when both the energy accumulation condition and the accelerated expansion condition are met, it is determined that the corresponding preliminary potential slip surface enters a critical unstable state, and it is output as a determined potential slip surface; if either of the above conditions is not met, it is determined that it has not entered a critical state.

[0113] The adjustment formula of the dynamic critical energy threshold value is as follows: ;

[0114] ;

[0115] ;

[0116] wherein: represents the kth slip surface dynamic critical energy threshold value at time t; represents the basic critical threshold value; represents the damage sensitivity coefficient; natural logarithm operator; represents the slip surface nearby time-varying amplitude of decrease; represents the initial integrity coefficient of the rock mass (reference value); represents the rainfall correction function; represents the cumulative rainfall; represents the rock mass integrity coefficient at the initial time ; represents the rock mass integrity coefficient at time t.

[0117] In this embodiment, the dynamic critical energy threshold value is reduced by increasing the structural damage correction term as the damage increases; the rainfall correction term is increased as the cumulative rainfall increases, which reduces the dynamic critical energy threshold value. At the same time, the area ratio and the expansion rate are added, and finally the high-risk area of continuous and accelerated expansion is determined to avoid isolated point interference.

[0118] Step S4. Obtain the strain energy density overage rate of the slip surface , the shear stress change rate with time , and the daily expansion rate index

[0119] , and the instability dynamics index is obtained by weighted summation based on the obtained indexes, and the specific expression is as follows: ;

[0120] wherein: represents the difference between the actual strain energy density of the slip surface and the dynamic critical energy threshold value; represents the maximum shear stress on the potential slip surface (real-time monitoring by the fiber Bragg grating stress sensing chain); represents the continuous area of the strain energy over-limit region; , and represent the weight coefficients; is the dynamic critical energy threshold value;

[0121] Different warning mechanisms are set based on the value, for example in the interval range of 0.3 to 0.6, yellow warning is performed; in the range of 0.6-0.8, orange warning is output; if greater than 0.8, a red early warning is output; different response measures are set for different early warning levels, and the setting of the response measures belongs to the conventional technical means in the art, and thus is not described in detail in the present application.

[0122] The above only describes the preferred embodiments of the present application and is not used to limit the present application, and any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A tunnel portal landslide prediction method based on joint inversion of shallow seismic and stress, characterized in that, Comprise the following steps: Step S1. Arranging a plurality of micro intelligent seismic source nodes in the target area to form a high-frequency seismic source array, and obtaining slope first arrival wave travel time data based on the high-frequency seismic source array, constructing a geological model based on the slope first arrival wave travel time data, and calculating rock mass integrity related parameters according to the elastic wave velocity distribution to form a rock mass integrity coefficient matrix; Step S2. Implanting a fiber grating stress sensing chain in the construction drilling to capture real-time change data of the main stress direction and the main stress level of the rock mass, taking the rock mass integrity coefficient matrix as prior knowledge, embedding it into the stress field solving process, establishing the correlation between seismic wave propagation and stress, and constructing the internal relationship between stress and wave velocity, and then coupling and solving to obtain slope stress field distribution data; Step S3. According to the stress data and strain data, the strain energy density distribution of each point on the potential sliding surface is calculated, whether the calculated strain energy density exceeds the predetermined threshold value and whether it is continuous and in pieces are determined, if it is continuous and in pieces and the expansion rate of the piece area exceeds the predetermined threshold value, it is determined as a potential sliding surface; Step S4. Obtain the strain energy density over-standard rate, the shear stress change rate with time and the daily expansion rate index of the sliding surface, and perform weighted summation based on the obtained indexes to obtain a instability dynamics index for risk assessment.

2. The tunnel portal landslide prediction method of the shallow seismic-stress joint inversion according to claim 1, characterized in that, Arranging an environmental sensor and a three-axis accelerometer in the plurality of micro intelligent seismic source nodes to collect rainfall intensity and construction vibration data; Adjusting the seismic source excitation frequency based on the collected rainfall intensity and construction vibration data, comprising: If the real-time rainfall intensity is greater than the predetermined threshold value, the seismic source excitation frequency is the product of the rainfall-frequency adjustment coefficient and the rainfall intensity plus the basic frequency; If the construction vibration intensity is greater than the predetermined threshold value, the seismic source excitation frequency is the product of the vibration-frequency adjustment coefficient and the vibration intensity plus the basic frequency; If neither of the above two cases, the seismic source excitation frequency is the set basic frequency.

3. The tunnel portal landslide prediction method of the shallow seismic-stress joint inversion according to claim 1, characterized in that, The formation of the rock mass integrity coefficient matrix comprises the following steps: Based on the AIC automatic picker, the slope first arrival wave travel time is extracted from the received signal, based on the extracted travel time, lithology parameters, a preset weight matrix and a regularization coefficient, a target function composed of a data fitting term and a model smoothing term is constructed, and then the three-dimensional grid P-wave velocity data is calculated by combining the anisotropy parameter correction slowness model through the fast marching method, and the residual error is updated by the LSQR algorithm and the anisotropy constraint is introduced when the residual error exceeds the limit in the iterative solving process until the residual error converges; Based on the obtained three-dimensional grid P-wave data, the complete rock mass experimental wave velocity and the rock mass structure classification standard, the square of the inversion wave velocity and the complete wave velocity ratio of each grid point is calculated to generate a rock mass integrity coefficient matrix; then the grid rock mass state 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 combined with the three-dimensional grid.

4. The tunnel portal landslide prediction method of claim 3, wherein, After the AIC automatic picker extracts the slope first arrival wave travel time, the extracted slope first arrival wave travel time needs to be corrected, comprising the following steps: The difference between the real-time temperature and the standard temperature, the difference between the real-time humidity and the standard humidity are calculated respectively, the difference between the real-time temperature multiplied by the temperature calibration coefficient is added to the difference between the humidity multiplied by the humidity calibration coefficient, and then 1 is added to obtain a correction factor, and the correction factor is multiplied by the slope body first arrival wave travel time to obtain the corrected slope body first arrival wave travel time. 5.The tunnel portal landslide prediction method of the shallow seismic-stress joint inversion according to claim 1, characterized in that, The internal relationship between the stress and the wave velocity is represented based on a constructed wave velocity-stress coupling equation group, wherein the wave velocity-stress coupling equation group comprises a motion balance equation, a constitutive equation constrained by an integrity coefficient matrix, and a geometric equation; The motion balance equation is used to describe the balance relationship between the rock mass density, the second-order time derivative of displacement, and the stress divergence and the volume force; The constitutive equation constrained by the integrity coefficient matrix is used to associate the stress and the strain, and the integrity coefficient matrix is introduced to make the elastic stiffness tensor dynamically change with the integrity coefficient matrix; The geometric equation defines the strain through the displacement gradient and establishes the geometric relationship between the displacement and the strain.

6. The tunnel portal landslide prediction method of claim 5, wherein, For the elastic stiffness tensor, the integrity coefficient is used for correction: The elastic stiffness tensor is obtained by multiplying the integrity coefficient by the rock mass stiffness and adding the value obtained by multiplying the damaged rock mass stiffness by 1 minus the integrity coefficient.

7. The tunnel portal landslide prediction method of claim 5, wherein the shallow seismic-stress joint inversion is performed by using a genetic algorithm. The slope stress field distribution data obtained by the coupling solution comprises the following steps: Based on the constructed wave velocity-stress coupling equation group, a regularization least square optimization algorithm is used to solve the objective function constructed by constructing a target function comprising a stress matching term and a monitoring data constraint term, until the objective function is minimized, to obtain a stable stress field distribution; The stress matching term is the square sum of the difference between the stress output by the finite element model and the measured value of the borehole gravity; and the monitoring data constraint term is the square sum of the difference between the inverted stress and the FBG monitoring stress, multiplied by a regularization weight. 8.The tunnel portal landslide prediction method of the shallow seismic-stress joint inversion according to claim 1, characterized in that, The step S3 comprises the following steps: Based on the structural surface basic parameters and the rock mass integrity coefficient matrix, a plurality of candidate slip surfaces are divided, and the normal shear stress ratio is calculated for each candidate slip surface; if the shear stress ratio is greater than or equal to a set threshold value, the corresponding candidate surface is taken as a preliminary potential slip surface; The components corresponding to the stress tensor and the strain tensor are multiplied two by two, and then all the product results are summed, and finally the sum result is multiplied by a predetermined coefficient to obtain the strain energy density of the corresponding spatial point at the current time, and based on this, all spatial points are traversed to generate three-dimensional strain energy density distribution data; The energy component along the gradient direction of the normal vector of the slip surface in the strain energy density is filtered through the Dirac function, and then the filtered energy component is integrated on the microelement area of the preliminary potential slip surface, the energy values of all microelement areas are accumulated, and the 2D shear strain energy density distribution of the preliminary potential slip surface is obtained. The area of the region where the shear strain energy density on each slip surface is greater than or equal to a dynamic critical energy threshold value is counted, and a ratio of the area to a total area of the slip surface is calculated; if the ratio is greater than or equal to a predetermined proportion threshold value, the energy accumulation condition is met; the rate of change of the area ratio with time is calculated, and if the rate of change is greater than a predetermined rate threshold value, the accelerated expansion condition is met; when both the energy accumulation condition and the accelerated expansion condition are met, it is determined that the corresponding preliminary potential slip surface enters a critical instability state, and the preliminary potential slip surface is output as a determined potential slip surface; if either of the conditions is not met, it is determined that the preliminary potential slip surface does not enter the critical state. 9.The tunnel portal landslide prediction method of the shallow seismic-stress joint inversion according to claim 8, characterized in that, The dynamic critical energy threshold value is adjusted based on the following steps: An adjustment function combining a structure damage correction term and a rainfall correction term is used to adjust the dynamic critical energy threshold value, wherein the structure damage correction term is obtained by taking the natural logarithm of the ratio of the time-varying amplitude of the rock mass integrity coefficient to the initial rock mass integrity coefficient, adding 1, and then multiplying by a damage sensitivity coefficient to obtain a threshold correction proportion caused by rock mass damage; The rainfall correction term adjusts the threshold value according to the cumulative rainfall amount through a hyperbolic tangent function; The dynamic critical energy threshold value of each potential slip surface at the current time is obtained by multiplying the basic threshold value by 1 plus the sum of the structure damage correction terms, and then multiplying by the rainfall correction term. 10.The tunnel portal landslide prediction method of the shallow seismic-stress joint inversion according to claim 1, characterized in that, The strain energy density overage rate is obtained by subtracting the strain energy critical threshold value from the actual strain energy density of the potential slip surface, and then dividing the difference by the strain energy critical threshold value; The strain energy critical threshold value is obtained based on the following steps: The deviation value of the current rock mass integrity from the regional average level is obtained by subtracting the regional rock mass integrity reference value from the real-time rock mass integrity coefficient; The threshold correction factor is obtained by multiplying the rock mass integrity deviation value by a correction coefficient and then adding 1, and the corrected strain energy critical threshold value is obtained by multiplying the reference strain energy threshold value by the threshold correction factor.

Images