Tunnel water inrush monitoring and early warning method under water-rock coupling effect

Abstract

The invention discloses a tunnel water inrush monitoring and early warning method under the water-rock coupling effect, and relates to the technical field of karst tunnels. Through the design of the method, the positions of karst caves and fractures of tunnels containing complex geological conditions such as karst caves and fractures under the water-rock coupling effect can be detected, the real-time change of water-containing fracture stress is detected, address disaster influence factors are analyzed based on machine learning to obtain evaluation indexes of the address disaster influence factors, and the evaluation indexes of the address disaster influence factors are analyzed. Therefore, a tunnel water inrush monitoring and early warning system is established, and a more accurate and rapid early warning method is provided for the tunnel water inrush problem.

Description

Technical Field

[0001] The present invention relates to the technical field of karst tunnels, and specifically to a method for monitoring and warning of water inrush in tunnels under the action of water-rock coupling. Background Technique

[0002] The focus of transportation infrastructure construction has shifted to the western mountainous areas and water-rich areas with extremely complex topographic and geological conditions. Therefore, a batch of tunnel (road) engineering projects with high water inrush risks, such as "large buried depth, high ground pressure, high water pressure, and large flow rate", will emerge, resulting in huge challenges for the safe construction of karst tunnels. Once a disaster occurs, it will bring significant casualties and property losses, and even lead to the suspension of the project.

[0003] Since the current research on karst tunnels is relatively weak, existing technologies mostly rely on single-type sensors (such as piezometers or displacement sensors), making it difficult to comprehensively capture the multi-factor coupling characteristics of water inrush disasters. Moreover, due to the extremely fast occurrence time of tunnel water inrush, traditional early warning relying on manual monitoring is slow to respond and has a long response time, unable to meet the suddenness requirements of water inrush disasters.

[0004] Therefore, a new solution to the above problems is needed. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for monitoring and warning of water inrush in tunnels under the action of water-rock coupling, and to propose a more accurate and rapid warning method for tunnel water inrush problems.

[0006] To achieve the above purpose, the present invention provides the following technical solution: A method for monitoring and warning of water inrush in tunnels under the action of water-rock coupling, at least including the following steps:

[0007] S100: Obtain the position cross-section diagram of tunnel karst caves and water-containing fissures based on the high-density resistivity method;

[0008] S200: Real-time sense the stress change of fractured rock masses;

[0009] S300: Prediction of rock mass water inrush risk based on machine learning, and set according to whether the prediction result reaches the trigger warning condition. When the trigger warning condition is reached, a warning is issued.

[0010] Further, the S1 at least includes the following steps:

[0011] S101: Collect tunnel rock mass data based on the high-density resistivity method;

[0012] S102: Data rejection and correction processing;

[0013] S103: Inversion calculation based on the least square iteration method;

[0014] S104: Infer the cross-sectional diagrams of karst caves and fissures in tunnel rock masses based on S103.

[0015] Further, the acquisition of tunnel rock mass data based on the high-density resistivity method means arranging electrodes at the position of the tunnel crown in the dipole-dipole mode and then collecting data. This method is sensitive to water-containing fissures, can clearly identify the change of resistivity gradient, and can enhance the abnormal signal by combining the polarization ratio parameter. The karst caves are distinguished as high resistance and the water-containing fissures as low resistance.

[0016] The principle of the dipole-dipole mode is as follows:

[0017] The lengths of the power supply dipole AB and the measurement dipole MN are both a, and the distance between them is n×a, where n is the isolation coefficient and takes positive integers.

[0018] The power supply and measurement dipoles are arranged on the same line, and the detection depth is positively correlated with n×a. The empirical formula is D = n×a×k, where k is the influence factor of the medium resistivity and takes 0.5 - 0.7 in the low-resistance area.

[0019] Measure the apparent resistivity ρ by calculating the potential difference through the distribution difference of the artificial current field in the underground medium. a , and the formula is:

[0020]

[0021] where K = 6πa is the device coefficient; ΔV is the potential difference between MN; I is the power supply current.

[0022] When arranging, the electrode spacing n×a is usually 1 - 5 meters, and the value range of the isolation coefficient n is 1 - 8.

[0023] Arrange a straight measuring line along the axis of the tunnel crown, number the electrodes at equal intervals (E1, E2,..., EN), and ensure that the grounding resistance < 5 kΩ.

[0024] Adopt longitudinal rolling measurement, that is, fix the AB dipole and sequentially collect the data of the measurement dipole MN at different intervals to form an inverted trapezoidal detection section.

[0025] Further, the data rejection and correction processing means preprocessing the original tunnel data collected by the high-density resistivity method, including rejecting bad data, data splicing, and distortion correction. S102 at least includes the following steps:

[0026] First is the rejection of bad data. Use the single-point relative error method to calculate the relative error of the data collected repeatedly at the same measuring point and screen out the abnormal values. The formula is:

[0027]

[0028] If δ i > δ阈值 The threshold value, usually taken as 5%-10%, is then determined as a bad point and eliminated; di is the i-th observation value; is the mean value of multiple observations;

[0029] Secondly, data stitching and correction are carried out. Using the non-linear deviation correction of the survey line, that is, when geometric distortion occurs to the survey line due to terrain undulation, a piecewise linear correction model is adopted:

[0030] x' = x + Δx(y)

[0031] y' = y + Δy(x)

[0032] By polynomial fitting of the coordinate deviations of homologous points in the overlapping area of adjacent survey lines, a correction mapping relationship is established to eliminate the geometric position deviation caused by terrain undulation;

[0033] Then the observed potential difference is converted into an equivalent value on a horizontal terrain:

[0034]

[0035] Map the data of the inclined terrain to a horizontal reference plane;

[0036] where x is the horizontal coordinate of the pixel point; y is the vertical coordinate of the pixel point; Δx(y) and Δy(x) are piecewise functions; ρ 真 is the theoretical resistivity under the horizontal terrain; ρ is regarded as the apparent resistivity under the actual terrain.

[0037] Finally, preprocessing of distortion correction is carried out. A two-dimensional equalization filter is designed by the least square method to suppress single-point and linear distortion:

[0038]

[0039] For the electric field distortion caused by poor electrode contact or non-uniform medium, a current density compensation model is established:

[0040]

[0041] By iteratively adjusting the k value, the electric field distribution after correction satisfies the continuity condition of the Poisson equation, thereby suppressing local electric field distortion and linear interference;

[0042] In the formula: is the image value after correction; represents that the range of weighted average is a rectangular window, and the window size is controlled by n and m, W ij is the filter coefficient, determined by minimizing the residual ε = ||d - Wd|| 2 d(x + j, y + i) represents the data or pixel value at the position (x + j, y + i) in the original image; k is the distortion coefficient; φ0 is the theoretical potential gradient; is the measured potential gradient.

[0043] Further, S103 is to perform inversion inference on the optimized data to infer the underground resistivity distribution, that is, using the current continuity equation and least squares optimization to complete the inversion calculation;

[0044] For forward calculation using the current continuity equation, assuming the medium resistivity is ρ(x, z) and the current source is I, the potential φ satisfies the Poisson equation:

[0045]

[0046] The Poisson equation is expressed as a linear equation system by the finite element method (FEM):

[0047] K(m)φ = q

[0048] In the formula: δ is the Dirac function; r s represents the position of the point current source; m is the model parameter, where the logarithmic transformation of the resistivity: m i = lnρ i ; K is the stiffness matrix, which depends on the model parameter; q is the source term vector;

[0049] The weighted least squares optimization inversion objective function is adopted, and the objective function is:

[0050]

[0051] In the formula, d obs is the observed data, that is, the logarithm of the potential difference or the apparent resistivity; F(m) is the forward response; W d is the data weighting matrix, usually a diagonal matrix, and the elements are the reciprocals of the data standard deviations; Wm is the model regularization matrix, which constrains the smoothness or structure of the model; λ is the regularization parameter, which balances data fitting and model complexity;

[0052] The detailed steps of the inversion inference are as follows:

[0053] First, the nonlinear problem is Taylor-expanded at the initial model m k and the first-order term is retained:

[0054] F(m k+1 ) ≈ F(m k ) + J k Δm k

[0055] In the formula: m k is the model of the current iteration; J k [[ID=**63]]is the Jacobian matrix, that is, the sensitivity matrix; Δm k is the update amount of the model; It should be noted that there seems to be a small error in the original text where "II" in the description of the current source should probably be "I". Also, in the translation, the numbering of some lines might be a bit off due to the format preservation requirements, but the overall translation content is correct. If you have any further questions or need more detailed explanations, please feel free to ask.

[0056] Substitute into the objective function to obtain the linearized optimization problem:

[0057] Δm k = argmin Δm ||W d (r k - J k Δm)|| 2 + λ||W m (m k + Δm)|| 2

[0058] where J k Δm is the predicted impact of the model update amount on the observed data; Δm is the model update amount;

[0059] Solve this equation to obtain the update amount of the model:

[0060]

[0061] In the formula: is the weighted matrix term of the Jacobian matrix, describing the sensitivity of the data to the model update amount; is the model regularization matrix term, used to control the model complexity and avoid overfitting; The product of the data residual and the Jacobian matrix represents the contribution of the data difference to the model update; is the model regularization term, restricting the model complexity;

[0062] r k = d obs - F(m k ) is the residual, that is, the difference between the current model mk and the observed data;

[0063] Secondly, it is the calculation of the Jacobian matrix. The element J ij of the Jacobian matrix J represents the sensitivity of the i-th data to the j-th model parameter: Obtain the spatial integrals of the forward field φ and the adjoint field ψ by calculation:

[0064]

[0065] In the formula: Ω is the model domain; is the derivative of the resistivity with respect to the model parameter; and are the gradients of the forward field and the adjoint field, representing the spatial variation of the field

[0066] Then it is the regularization design. The regularization matrix W m usually adopts a smoothing constraint:

[0067]

[0068] wherein and are the second-order derivatives of the model in the x and y directions (i.e., the Laplacian operator), which represent the spatial smoothness of the model

[0069] Finally, set the iteration termination conditions, which are divided into residual convergence and model change convergence;

[0070] Residual convergence is that the relative root mean square error (RMS) is less than the threshold:

[0071]

[0072] where: d obs,i is the i-th observed data; F i (m) is the forward response of the i-th observed data, that is, the observed value predicted by the model m; σ i is the standard deviation of the i-th observed data; N is the total number of observed data;

[0073] Model change convergence is ||Δm k || < ∈.

[0074] Furthermore, the real-time perception of the stress change of the fractured rock mass is carried out by a real-time perception device for the stress change of the fractured rock mass, and the real-time perception device for the stress change of the fractured rock mass includes a stress sensor, a stress transmission line, a host, a rock mass, a water-containing fracture and a tunnel;

[0075] The application of the real-time perception device for the stress change of the fractured rock mass at least includes the following steps:

[0076] First, according to the position of the cross-sectional diagram of the water-containing fracture obtained by S100, insert the stress sensor borehole into the rock mass;

[0077] The depth inserted into the rock mass should reach the water-containing fracture to measure the stress σ n and water pressure P of the rock mass;

[0078] Transmit the measured stress data to the host through the stress transmission line. Secondly, set the ultimate stress value [σ u and the ultimate shear strength of the fractured rock mass in the host

[0079] Finally, according to the real-time stress change of the monitored rock mass, analyze the effective stress σ' n and the effective shear stress τ' n of the water-containing fracture by the water-rock coupling method, and judge the degree of fracture cracking in combination with fracture mechanics;

[0080] The ultimate stress value [σ u and the ultimate shear strength of the fractured rock mass It is to take on-site samples of the tunnel, use a uniaxial test instrument to test the uniaxial compressive strength of the rock mass, so as to obtain the compressive strength of the fractured rock mass, and use this as the ultimate stress value [σ u and the ultimate shear strength Then the obtained values are input into the host through a mechanical key panel.

[0081] Furthermore, the rock mass water inrush risk prediction based on machine learning at least includes the following steps:

[0082] First, collect the geological disaster classification data of the rock mass water inrush caused by the tunnel. Through the stress sensor, collect the changes in the tunnel rock mass stress σ(t) and seepage pressure P(t) in the past half month, collect the crown settlement D(t) at the top of the tunnel in the past half month, collect the historical maximum rainfall in the past 10 years according to the information of the meteorological bureau, conduct data collection, and further transmit the collected data to the data preprocessing module;

[0083] Secondly, organize, screen and optimize the data. In the data preprocessing module, use the Z-score method to detect outliers for the rock mass stress σ(t), seepage pressure P(t), and crown settlement D(t). Its expression is:

[0084]

[0085] In the formula: Z i is the Z-score of the i-th data point; x i is the value of the i-th data point; μ is the mean of the data set; e is the standard deviation of the data set;

[0086] Normalize the data at the outlet. The expression is:

[0087]

[0088] In the formula: x is the original data collected; min(x) is the minimum value in the data set; max(x) is the maximum value in the data set;

[0089] For the rainfall feature extraction, extract the statistical features of the historical maximum rainfall R max . First, calculate the ratio of the current period rainfall R(t) to the historical maximum value

[0090]

[0091] Then introduce the attenuation factor λ to weight the influence of historical rainfall:

[0092]

[0093] Among them, R(t-k) is the past rainfall with a time step of k;

[0094] Multi - factor weight distribution, construct a data matrix for the index rock mass stress, osmotic pressure, settlement, rainfall ratio, and cumulative rainfall

[0095]

[0096] Among them, x ij represents the value of the j - th index at the i - th time point, n is the number of time steps, and m is the number of indices, indicating how many different indices there are;

[0097] Calculate the information entropy of the j - th index:

[0098]

[0099] Among them, is the normalization constant to ensure that the value of the information entropy is between [0, 1]; is the normalized probability, representing the relative proportion of the j - th index at the i - th time point; P ij is the standardized value of the j - th index at the i - th time point; is the sum of the j - th index over all time steps; E j is the information entropy, measuring the uncertainty of the j - th index;

[0100] Then calculate the index weight W j , reflecting its contribution:

[0101]

[0102] In the formula: m is the number of indices;

[0103] Immediately construct a risk probability model, define the water - inrush event label Y ∈ {0, 1}, and the input feature vector is X = [σ norm , P norm , D norm , R ratio , R cum , σ norm is the standardized value of the rock mass stress; P norm is the standardized value of the osmotic pressure; D norm is the standardized value of the settlement; R ratio is the rainfall ratio; R cum is the cumulative rainfall;

[0104] The logistic regression model is: [[ID=6৮]]

[0105]

[0106] In the formula: Y is the predicted target variable, indicating whether a water inrush disaster occurs. Y = 1 means it occurs, and Y = 0 means it does not occur; β0 is the intercept term; βj is the feature coefficient; is a linear combination in the logistic regression model, which is transformed into probability through an exponential function;

[0107] Finally, a risk probability model is constructed for predicting the risk of water inrush in rock masses. The log-likelihood function is constructed by the maximum likelihood method and the parameters are optimized and solved:

[0108]

[0109] where, P i = P(Y = 1|x i ) is the predicted probability of the i-th sample, that is, the water inrush probability calculated through the logistic regression model; y i is the actual label of the i-th sample; the parameter estimation value that maximizes lnL(β) is solved by iterative method of Newton-Raphson

[0110] When the risk model is constructed, the early warning dynamic threshold is set accordingly. First, the threshold is calibrated. Based on the historical water inrush event data, the optimal threshold θ is determined through the ROC curve, and it needs to satisfy:

[0111] θ = argmax(Youden index) = argmax(sensitivity + specificity - 1)

[0112] Set the early warning dynamic threshold and calculate the risk probability P new = P(Y = 1|x new ) in real time, and trigger the early warning condition:

[0113]

[0114] Compared with the prior art, the beneficial effects of the present invention are:

[0115] Through the design of the method, the present invention can detect the positions of karsts and fissures in a tunnel with complex geological conditions such as karsts and fissures under the action of water-rock coupling, detect the real-time change of the stress of water-containing fissures, analyze the influencing factors of geological disasters based on machine learning to obtain their evaluation indexes, and thus establish a tunnel water inrush monitoring and early warning system, providing a more accurate and rapid early warning method for the tunnel water inrush problem. BRIEF DESCRIPTION OF THE DRAWINGS

[0116] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for the description of the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention, and those of ordinary skill in the art can also obtain other drawings based on these drawings without creative efforts.

[0117] Figure 1 is the flow chart of the present invention;

[0118] Figure 2 is the flow chart of the present invention for obtaining the position cross-section diagram of tunnel karst caves and water-bearing fissures based on the high-density electrical method;

[0119] Figure 3 is the schematic diagram of the device for real-time sensing of stress changes in fractured rock masses of the present invention;

[0120] Figure 4 is the stress analysis diagram of the water-bearing fissure of the present invention.

[0121] In the figure: 1, stress sensor; 2, stress transmission line; 3, host; 4, rock mass; 5, water-bearing fissure; 6, tunnel. Specific embodiments

[0122] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments.

[0123] Please refer to Figure 1 , the method for monitoring and warning of tunnel water inrush under the action of water-rock coupling includes at least the following steps:

[0124] S100: Obtain the position cross-section diagram of tunnel karst caves and water-bearing fissures based on the high-density electrical method;

[0125] S200: Real-time sense the stress change of fractured rock masses;

[0126] S300: Prediction of rock mass water inrush risk based on machine learning, which is set according to whether the prediction result reaches the trigger warning condition, and gives a warning when the trigger warning condition is reached.

[0127] In this embodiment, an additional step can be added to obtain a tunnel water inrush monitoring and warning system;

[0128] S400: Tunnel water inrush monitoring and warning system under the action of multi-factor coupling.

[0129] S1 includes at least the following steps:

[0130] S101: Collect tunnel rock mass data based on the high-density electrical method;

[0131] S102: Data rejection and correction processing;

[0132] S103: Inverse calculation based on the least square iteration method;

[0133] S104: Infer the cross-sectional diagram of karst caves and fissures in tunnel rock mass based on S103.

[0134] Collect tunnel rock mass data based on the high-density resistivity method, that is, arrange electrodes at the position of the tunnel crown in the dipole-dipole mode and then collect data. This method is sensitive to water-containing fissures, can clearly identify the change of resistivity gradient, and can enhance the abnormal signal by combining the polarization ratio parameter. Distinguish according to high resistance for karst caves and low resistance for water-containing fissures.

[0135] The principle of the dipole-dipole mode is:

[0136] The lengths of the power supply dipole AB and the measurement dipole MN are both a, and the distance between them is n×a, where n is the isolation coefficient and takes positive integers.

[0137] The power supply and measurement dipoles are arranged on the same line, and the detection depth is positively correlated with n×a. The empirical formula is D = n×a×k, where k is the influence factor of the medium resistivity, and takes 0.5 - 0.7 in the low-resistance area.

[0138] Measure the potential difference to calculate the apparent resistivity ρ through the distribution difference of the artificial current field in the underground medium a , and the formula is:

[0139]

[0140] Among them, K = 6πa is the device coefficient; ΔV is the potential difference between MN; I is the power supply current;

[0141] When arranging, the electrode spacing n×a is usually 1 - 5 meters, and the value range of the isolation coefficient n is 1 - 8;

[0142] Arrange a straight survey line along the axis of the tunnel crown, and number the electrodes at equal intervals (E1, E2,..., EN) to ensure that the grounding resistance < 5kΩ;

[0143] Adopt longitudinal rolling measurement, that is, fix the AB dipole and collect the data of the measurement dipole MN at different intervals in turn to form an inverted trapezoidal detection section.

[0144] Data rejection and correction processing means preprocessing the original tunnel data collected by the high-density resistivity method, including rejecting bad data, data splicing and distortion correction. S102 at least includes the following steps:

[0145] First is bad data rejection. Use the single-point relative error method to calculate the relative error of the data collected repeatedly at the same measuring point, and screen out the abnormal values. The formula is:

[0146]

[0147] If δ i > δ 阈值The threshold value, usually taken as 5%-10%, is then determined as a bad pixel and eliminated; di is the ith observed value; is the mean value of multiple observations;

[0148] Secondly is data stitching and calibration. Using the non-linear deviation calibration of the survey line, that is, when geometric distortion occurs due to terrain undulation of the survey line, a piecewise linear calibration model is adopted:

[0149] x' = x + Δx(y)

[0150] y' = y + Δy(x)

[0151] By polynomial fitting of the coordinate deviation of the same-name points in the overlapping area of adjacent survey lines, a calibration mapping relationship is established to eliminate the geometric position deviation caused by terrain undulation;

[0152] Then the observed potential difference is converted into an equivalent value on a horizontal terrain:

[0153]

[0154] Map the data of the inclined terrain to the horizontal reference plane;

[0155] Among them, x is the horizontal coordinate of the pixel point; y is the vertical coordinate of the pixel point; Δx(y) and Δy(x) are piecewise functions; ρ 真 is the theoretical resistivity under the horizontal terrain; ρ is regarded as the apparent resistivity under the actual terrain.

[0156] Finally is the preprocessing of distortion correction. A two-dimensional equalization filter is designed using the least squares method to suppress single-point and linear distortions:

[0157]

[0158] For the electric field distortion caused by poor electrode contact or non-uniform medium, a current density compensation model is established:

[0159]

[0160] By iteratively adjusting the k value, the electric field distribution after calibration satisfies the continuity condition of the Poisson equation, thereby suppressing local electric field distortion and linear interference;

[0161] In the formula: is the image value after calibration; represents that the range of weighted average is a rectangular window, and the window size is controlled by n and m, W ij is the filter coefficient, by minimizing the residual ε = ||d - Wd|| 2 determined; d(x + j, y + i) represents the data or pixel value at the position (x + j, y + i) in the original image; k is the distortion coefficient; φ0 is the theoretical potential gradient; is the measured potential gradient.

[0162] S103 is to perform inversion inference on the optimized data to infer the underground resistivity distribution, that is, using the current continuity equation and least squares optimization to complete the inversion calculation;

[0163] Using the current continuity equation for forward calculation, assuming the medium resistivity is ρ(x,z) and the current source is II, then the potential φ satisfies the Poisson equation:

[0164]

[0165] The Poisson equation is expressed as a linear equation system by the finite element method (FEM):

[0166] K(m)φ = q

[0167] In the formula: δ is the Dirac function; r s represents the position of the point current source; m is the model parameter, where the logarithmic transformation of the resistivity: m i = lnρ i ; K is the stiffness matrix, which depends on the model parameter; q is the source term vector;

[0168] Adopt the weighted least squares optimization inversion objective function, and the objective function is:

[0169]

[0170] In the formula, d obs is the observed data, that is, the logarithm of the potential difference or the apparent resistivity; F(m) is the forward response; W d is the data weighting matrix, usually a diagonal matrix, and the elements are the reciprocals of the data standard deviations; Wm is the model regularization matrix, which constrains the smoothness or structure of the model; λ is the regularization parameter, which balances data fitting and model complexity;

[0171] The detailed steps of the inversion inference are:

[0172] First, expand the nonlinear problem at the initial model m k and retain the first-order term:

[0173] F(m k+1 ) ≈ F(m k ) + J k Δm k

[0174] In the formula: m k is the model of the current iteration; J k is the Jacobian matrix, that is, the sensitivity matrix; Δm k is the update amount of the model;

[0175] Substitute into the objective function to obtain the linearized optimization problem:

[0176] Δm k =argmin Δm ||W d (r k -J k Δm)|| 2 +λ||W m (m k +Δm)|| 2

[0177] where J k Δm is the predicted impact of the model update amount on the observed data; Δm is the model update amount;

[0178] Solve this equation to obtain the update amount of the model:

[0179]

[0180] In the formula: is the weighted matrix term of the Jacobian matrix, describing the sensitivity of the data to the model update amount; is the model regularization matrix term, used to control the model complexity and avoid overfitting; The product of the data residual and the Jacobian matrix represents the contribution of the data difference to the model update; is the model regularization term, restricting the model complexity;

[0181] r k =d obs -F(m k ) is the residual, that is, the difference between the current model mk and the observed data;

[0182] Secondly, it is the calculation of the Jacobian matrix. The element J ij of the Jacobian matrix J represents the sensitivity of the i-th data to the j-th model parameter: Obtain the spatial integrals of the two fields by calculating the forward field φ and the adjoint field ψ:

[0183]

[0184] In the formula: Ω is the model domain; is the derivative of the resistivity with respect to the model parameter; and are the gradients of the forward field and the adjoint field, representing the spatial variation of the field

[0185] Then it is the regularization design. The regularization matrix W m usually adopts a smoothing constraint:

[0186]

[0187] in and is the second-order derivative of the model in the x and y directions (i.e., the Laplace operator), which represents the spatial smoothness of the model

[0188] Finally, set the iteration termination conditions, which are divided into residual convergence and model change convergence;

[0189] The residual converges to the relative root mean square error (RMS) less than the threshold:

[0190]

[0191] Where: d obs,i is the i-th observation data; F i (m) is the forward response of the i-th observation data, that is, the observation value predicted by model m; σ i is the standard deviation of the i-th observation data; N is the total number of observation data;

[0192] The model change converges to ||Δm k ||<∈.

[0193] Real-time sensing of stress changes in fractured rock mass is performed by a real-time sensing device for stress changes in fractured rock mass, which includes a stress sensor 1, a stress transmission line 2, a host 3, a rock mass 4, a water-bearing fracture 5, and a tunnel 6;

[0194] The application of the device for real-time sensing of stress changes in fractured rock mass includes at least the following steps:

[0195] First, according to the position of the water-bearing fracture cross-section obtained in S100, the stress sensor 1 is drilled and inserted into the rock mass (4);

[0196] The depth of the rock mass 4 is to reach the water-bearing fissure 5 and measure the stress σ of the rock mass. n and water pressure P;

[0197] The measured stress data is transmitted to the host computer 3 through the stress transmission line 2, and then the ultimate stress value of the fractured rock mass [σ u ] and ultimate shear strength

[0198] Finally, according to the real-time changes of rock mass stress, the effective stress σ' of water-bearing fractures is analyzed by water-rock coupling method. n and effective shear stress τ' n , combined with fracture mechanics to determine the degree of cracking;

[0199] Ultimate stress value of fractured rock mass [σ u ] and ultimate shear strength Take on-site samples of Tunnel 6, and use a uniaxial test instrument to measure the uniaxial compressive strength of the rock mass, so as to obtain the compressive strength of the fractured rock mass, and use this as the ultimate stress value [σ u and the ultimate shear strength Then input the obtained values into the host 3 through the mechanical key panel;

[0200] The effective stress σ' of the water-bearing fracture is analyzed by the water-rock coupling method n and the effective shear stress τ' n , as Figure 4 shown, the effective stress of the rock mass fracture is obtained according to the effective stress law of the saturated porous elastic medium: σ' n =σ n -αP;

[0201] In the formula: σ n is the total stress of the rock mass, measured by the stress sensor; P is the water pressure of the water-bearing fractured rock mass, also measured by the stress sensor, and its calculation formula is: P = r w (H - z), r w is the unit weight of water, H is the water head, z is its position elevation; α is the effective stress coefficient of the fracture.

[0202] For a completely saturated rough rock fracture, when σ n ≥p, When σ n <p, α = 1;

[0203] In the formula: k ao , u nmax are the initial stiffness and the maximum normal closure of the fracture respectively.

[0204] The influence of the water-rock coupling action on the stability of the water-bearing fracture of the rock mass, according to the above expression of the effective stress, the classical Coulomb criterion is modified as follows: τ' n =c t +μ t (σ n -αP);

[0205] In the formula: c t , μ t are the friction factor and cohesion of a single fracture considering the time effect of the water-rock interaction respectively.

[0206] Judge whether the crack is cracked according to fracture mechanics, and the specific steps are as follows:

[0207] The crack propagation belongs to the compression-reducing crack opening, which belongs to the pure type II crack problem. Then, considering the water-rock coupling action, the calculation formula of the stress intensity factor at the tip of the water-bearing fracture of the tunnel surrounding rock is: In the formula, a is the crack length;

[0208] When the rock mass stress intensity factor reaches the crack initiation point, K Ⅱ >K ⅡC This is used to determine whether the crack is open.

[0209] The prediction of rock mass water inrush risk based on machine learning includes at least the following steps:

[0210] First, we collect classified data on geological disasters caused by rock inrush in tunnels. We use stress sensor 1 to collect the changes in tunnel rock stress σ(t) and seepage pressure P(t) over the past two weeks. We also collect the vault settlement D(t) over the past two weeks at the top of the tunnel. We also collect the historical maximum rainfall over the past 10 years based on information from the Meteorological Bureau. The collected data is then transmitted to the data preprocessing module.

[0211] Secondly, the data is sorted, screened and optimized. In the data preprocessing module, the rock stress σ(t), seepage pressure P(t) and vault settlement D(t) are detected using the Z-score method. The expression is:

[0212]

[0213] Where: Z i is the Z-score of the i-th data point; x i is the value of the i-th data point; μ is the mean of the data set; e is the standard deviation of the data set;

[0214] Normalize the data and export it. The expression is:

[0215]

[0216] Where: x is the collected raw data; min(x) is the minimum value in the data set; max(x) is the maximum value in the data set;

[0217] For rainfall feature extraction, extract the historical maximum rainfall R max Statistical characteristics of the first calculation of the ratio of the current period rainfall R (t) to the historical maximum value

[0218]

[0219] Then the attenuation factor λ is introduced to weight the impact of historical rainfall:

[0220]

[0221] Where R(tk) is the past rainfall with time step k;

[0222] Multi-factor weight allocation, constructing a data matrix for the indicators of rock mass stress, osmotic pressure, settlement, rainfall ratio, and cumulative rainfall

[0223]

[0224] Among them, x ij represents the value of the j-th indicator at the i-th time point, n is the number of time steps, and m is the number of indicators, indicating how many different indicators there are;

[0225] Calculate the information entropy of the j-th indicator:

[0226]

[0227] Among them, is the normalization constant to ensure that the value of the information entropy is between [0,1]; is the normalized probability, representing the relative proportion of the j-th indicator at the i-th time point; P ij is the standardized value of the j-th indicator at the i-th time point; is the sum of the j-th indicator over all time steps; E j is the information entropy, measuring the uncertainty of the j-th indicator;

[0228] Then calculate the indicator weight W j , reflecting its contribution:

[0229]

[0230] In the formula: m is the number of indicators;

[0231] Immediately construct a risk probability model, define the water inrush event label Y ∈ {0,1}, and the input feature vector is X = [σ norm , P norm , D norm , R ratio , R cum , σ norm is the standardized value of rock mass stress; P norm is the standardized value of osmotic pressure; D norm is the standardized value of settlement; R ratio is the rainfall ratio; R cum is the cumulative rainfall;

[0232] The logistic regression model is:

[0233]

[0234] In the formula: Y is the predicted target variable, indicating whether a water inrush disaster occurs, Y = 1 means it occurs, and Y = 0 means it does not occur; β0 is the intercept term; βj is the feature coefficient; It is a linear combination in the logistic regression model and is transformed into a probability through an exponential function;

[0235] Finally, a risk probability model is constructed for predicting the water inrush risk of rock masses. The log-likelihood function is constructed by the maximum likelihood method and the parameters are optimized and solved:

[0236]

[0237] Among them, P i =P(Y = 1|x i ) is the predicted probability of the i-th sample, that is, the water inrush probability calculated by the logistic regression model; y i is the actual label of the i-th sample; the parameter estimate value that maximizes lnL(β) is solved by iteration using the Newton-Raphson method

[0238] After the risk model is constructed, the early warning dynamic threshold is set accordingly. First, the threshold is calibrated. Based on the historical water inrush event data, the optimal threshold θ is determined through the ROC curve and needs to satisfy:

[0239] θ = argmax(Youden index) = argmax(sensitivity + specificity - 1)

[0240] Set the early warning dynamic threshold and calculate the risk probability P new =P(Y = 1|x new ) in real time, and trigger the early warning condition:

[0241]

[0242] Multi-factor coupling action tunnel water inrush monitoring and early warning system:

[0243] The early warning of water inrush from water-bearing fissures under the action of water-rock interaction is combined with the early warning of the rock mass water inrush risk prediction model based on machine learning to form a tunnel water inrush early warning system;

[0244] Establishment of early warning for the expansion of water-bearing fissures under the action of water-rock interaction;

[0245] First, according to the rock mass stress value transmitted by the stress sensor 1 to the host 3, and then by analyzing the effective normal stress σ' n and shear stress τ' n of the water-bearing fissures under the action of water-rock interaction, it is judged whether the fissures expand according to fracture mechanics, and the water inrush early warning is formed based on the cracking degree. Secondly, it is judged whether a single fissure expands according to fracture mechanics, and finally the water inrush early warning is formed based on the cracking degree;

[0246] Setting of water inrush early warning threshold: Level 1 early warning (yellow): K Ⅱ ≥0.8K ⅡC, indicating that the crack enters the subcritical propagation stage; Secondary warning (orange): K Ⅱ ≥0.9K ⅡC , emergency reinforcement measures need to be initiated; Tertiary warning (red): K Ⅱ ≥0.9K ⅡC , immediately evacuate personnel and close the working face;

[0247] Closely combine the water - inrush warning of water - containing fissure expansion under water - rock interaction with the warning of the rock mass water - inrush risk prediction model based on machine learning to form a tunnel water - inrush monitoring and warning system under multi - coupling action:

[0248] The stress change of the rock mass detected by the stress sensor 1 is transmitted to the host 3. In the host, the effective stress change of the fractured surrounding rock under water - rock coupling action is analyzed, and the ultimate shear strength of the rock mass obtained through uniaxial test and the stress allowable value [σ u are input into the host 3 through the mechanical keyboard. The fracture mechanics is used to judge the degree of fissure cracking. The set warning condition threshold is also input into the host 3 through the mechanical keyboard. At the same time, the tunnel water - inrush monitoring and warning threshold under multi - factor coupling action is also input into the host 3 through the mechanical keyboard. Then, alarm sound and light devices are evenly arranged on both sides of the tunnel arch waist, and the alarm sound and light devices are connected to the host 3;

[0249] When a primary warning occurs, that is, a potential risk, it is pushed to the monitoring terminal through the API to generate a yellow alarm log; When a secondary warning occurs, that is, a high - probability fault, the audible and visual alarm device is triggered, and a text message is sent to the responsible person synchronously; For the tertiary warning, that is, an emergency fault, the emergency shutdown protocol is automatically started, the alarm sound and light device is turned on for response, and the fault snapshot data is uploaded to the cloud.

[0250] It should be noted that the same computer terminal is connected to the high - density detection of the position of karst cave water and rock mass water - containing fissures, the real - time perception of the stress change of fractured rock mass, the rock mass water - inrush risk prediction based on machine learning, and the tunnel water - inrush monitoring and warning system under factor coupling action. First, the data obtained based on the high - density electrical method needs to be transmitted to the computer terminal for processing, and the cross - section map of the tunnel fissure position is obtained after processing; Secondly, the stress change of the fractured rock mass is perceived in real - time. According to the position of the water - containing fissure obtained from the cross - section map, the stress sensor is inserted into the position of the fissure, and the obtained rock mass stress data is input into the host 3. Similarly, the host 3 is transmitted to the computer terminal through data; Finally, for the rock mass water - inrush risk prediction based on machine learning and the tunnel water - inrush monitoring and warning system under factor coupling action, after constructing each risk model in the computer terminal, the warning dynamic threshold is set accordingly, which also includes the stress change of the rock mass fissure. First, the threshold is calibrated, and then the corresponding warning is set according to each factor. The computer terminal controls the alarm sound and light device to give a warning according to the tunnel danger situation.

[0251] For those skilled in the art, it is obvious that the present invention is not limited to the details of the above-described exemplary embodiments, and the present invention can be implemented in other specific forms without departing from the spirit or essential features of the present invention. Therefore, in any respect, the embodiments should be regarded as exemplary and non-limiting. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, all changes falling within the meaning and scope of the equivalent elements of the claims are intended to be embraced within the present invention. Any reference signs in the claims should not be construed as limiting the claims involved.

Claims

1. A method for monitoring and warning of water inrush in tunnels under the action of water-rock coupling, characterized in that: It includes at least the following steps: S100: Obtain the position cross-section diagrams of tunnel karsts and water-bearing fissures based on the high-density electrical method; S200: Real-time sense the stress changes of fractured rock masses; S300: Prediction of the risk of water inrush from rock masses based on machine learning, which is set according to whether the prediction result reaches the trigger warning condition, and give a warning when the trigger warning condition is reached.

2. The method for monitoring and warning of tunnel water inrush under the action of water-rock coupling according to claim 1, characterized in that: The S1 at least includes the following steps: S101: Collect tunnel rock mass data based on the high-density electrical method; S102: Data rejection and correction processing; S103: Conduct inversion calculation based on the least squares iteration method; S104: Infer the cross-section diagrams of karsts and fissures in the tunnel rock mass based on S103.

3. The method for monitoring and warning of water inrush in tunnels under the action of water-rock coupling according to claim 2, characterized in that: Collecting tunnel rock mass data based on the high-density electrical method means arranging electrodes at the position of the tunnel crown in the dipole-dipole manner and then collecting data. This method is sensitive to water-bearing fissures, can clearly identify the resistivity gradient changes, and combines the polarization ratio parameter to enhance the abnormal signal, and distinguishes according to high resistance for karsts and low resistance for water-bearing fissures; The principle of the dipole-dipole method is: The lengths of the power supply dipole AB and the measurement dipole MN are both a, and the distance between them is n×a, where n is the isolation coefficient and takes a positive integer; The power supply and measurement dipoles are arranged on the same line, and the detection depth is positively correlated with n×a. The empirical formula is D = n×a×k, where k is the medium resistivity influence factor, and 0.5 - 0.7 is taken for the low-resistance area; Measure the potential difference and calculate the apparent resistivity ρ through the distribution difference of the artificial current field in the underground medium a , and the formula is as follows: Among them, K = 6πa is the device coefficient; ΔV is the potential difference between MN; I is the power supply current; Arrange a straight survey line along the axis of the tunnel crown, number the electrodes at equal intervals, and ensure that the grounding resistance <5 kΩ; Adopt longitudinal rolling measurement, that is, fix the AB dipole and sequentially collect the data of the measurement dipole MN at different intervals to form an inverted trapezoidal detection cross-section.

4. The tunnel water inrush monitoring and early warning method under the action of water-rock coupling according to claim 3, characterized in that: The data rejection and correction processing is to perform preprocessing of rejecting bad data, data splicing, and distortion correction on the original tunnel data collected by the high-density electrical method. The S102 at least includes the following steps: First is the rejection of bad data. Calculate the relative error of the data collected repeatedly at the same measurement point by using the single-point relative error method, and screen out the abnormal values. The formula is: If δ i > δ 阈值 threshold, usually taken as 5%-10%, it is determined as a bad point and eliminated; di is the i-th observed value; is the mean of multiple observations; Secondly is data splicing and correction. Use the non-linear deviation correction of the survey line. That is, when geometric distortion occurs to the survey line due to terrain undulation, adopt the piecewise linear correction model: x' = x + Δx(y) y' = y + Δy(x) Establish a correction mapping relationship by polynomial fitting of the coordinate deviations of the same-name points in the overlapping area of adjacent survey lines to eliminate the geometric position deviation caused by terrain undulation; Then convert the observed potential difference into the equivalent value of the horizontal terrain: Map the inclined terrain data to the horizontal reference plane; where x is the horizontal coordinate of the pixel point; y is the vertical coordinate of the pixel point; Δx(y) and Δy(x) are piecewise functions; ρ 真 is the theoretical resistivity under horizontal terrain; ρ is regarded as the apparent resistivity under actual terrain; Finally is the preprocessing of distortion correction. Design a two-dimensional equalization filter by using the least squares method to suppress single-point and linear distortions: For the electric field distortion caused by poor electrode contact or non-uniform medium, establish a current density compensation model: Iteratively adjust the k value to make the corrected electric field distribution satisfy the continuity condition of the Poisson equation, so as to suppress local electric field distortion and linear interference; Wherein: is the corrected image value; indicates that the range of weighted average is a rectangular window, and the window size is controlled by n and m, and W ij is the filter coefficient, which is determined by minimizing the residual ε = ||d - Wd|| 2 ; d(x + j, y + i) represents the data or pixel value at the position (x + j, y + i) in the original image; k is the distortion coefficient; φ0 is the theoretical potential gradient; and ▽φ is the measured potential gradient.

5. The method for monitoring and warning of water inrush in tunnels under the action of water-rock coupling according to claim 4, characterized in that: The S103 is to perform inversion inference on the optimized data to infer the underground resistivity distribution, that is, use the current continuity equation and least squares optimization to complete the inversion calculation; Forward calculation is carried out using the current continuity equation. Assuming that the resistivity of the medium is ρ(x, z) and the current source is I, the potential φ satisfies the Poisson equation: The Poisson equation is expressed as a linear equation system by the finite element method: K(m)φ = q where: δ is the Dirac function; r s represents the position of the point current source; m is the model parameter, where the logarithmic transformation of the resistivity: m i = lnρ i ; K is the stiffness matrix, which depends on the model parameter; q is the source term vector; The weighted least squares optimization inversion objective function is adopted, and the objective function is: where d obs is the observed data, i.e., the logarithm of the potential difference or the apparent resistivity; F(m) is the forward response; W d is the data weighting matrix, usually a diagonal matrix with elements being the reciprocals of the data standard deviations; Wm is the model regularization matrix that constrains the smoothness or structure of the model; λ is the regularization parameter that balances data fitting and model complexity; The detailed steps of the inversion inference are as follows: First, the non - linear problem is Taylor expanded at the initial model m k and the first - order term is retained: F(m k+1 )≈F(m k )+J k Δm k where: m k is the model of the current iteration; J k is the Jacobian matrix, i.e., the sensitivity matrix; Δm k is the update amount of the model; Substitute into the objective function to obtain the linearized optimization problem: Δm k = argmin Δm ||W d (r k - J k Δm)|| 2 + λ||W m (m k + Δm)|| 2 Among them, J k Δm is the predicted influence of the model update amount on the observed data; Δm is the model update amount; Solve this equation to obtain the update amount of the model: In the formula: is the weighted matrix term of the Jacobian matrix, which describes the sensitivity of the data to the model update amount; is the model regularization matrix term, which is used to control the model complexity and avoid overfitting; The product of the data residual and the Jacobian matrix, which represents the contribution of the data difference to the model update; is the model regularization term, which restricts the model complexity; r k = d obs - F(m k ) is the residual, that is, the difference between the current model mk and the observed data; Next is the calculation of the Jacobian matrix. The element J of the Jacobian matrix J ij represents the sensitivity of the i-th data to the j-th model parameter: The spatial integrals of the forward field φ and the adjoint field ψ are obtained by calculation: where: Ω is the model domain; is the derivative of resistivity with respect to the model parameters; ▽φ and ▽ψ are the gradients of the forward field and the adjoint field, representing the spatial variation of the field Then comes the regularization design, the regularization matrix W m Smoothing constraints are usually adopted: where and are the second-order derivatives of the model in the x and y directions (i.e., the Laplacian operator), which represent the spatial smoothness of the model Finally, set the iteration termination conditions, which are divided into residual convergence and model change convergence; Residual convergence means that the relative root mean square error is less than the threshold: where: d obs,i is the i-th observed data; F i (m) is the forward response of the i-th observed data, i.e., the observed value predicted by model m; σ i is the standard deviation of the i-th observed data; N is the total number of observed data; The model change converges to ||Δm k || < ∈.

6. The monitoring and early warning method for tunnel water inrush under the action of water-rock coupling according to claim 1, characterized in that: The real-time perception of the stress change of fractured rock mass is carried out by a real-time perception device for the stress change of fractured rock mass. The real-time perception device for the stress change of fractured rock mass includes a stress sensor (1), a stress transmission line (2), a host (3), a rock mass (4), a water-containing fracture (5) and a tunnel (6); The application of the real-time perception device for the stress change of fractured rock mass at least includes the following steps: First, according to the position of the water-containing fracture cross-section diagram obtained by S100, insert the stress sensor (1) into the rock mass (4) by drilling; The depth of incorporation into the rock mass (4) should reach the water-bearing fracture (5) to measure the stress σ of the rock mass n and the water pressure P; The measured stress data is transmitted into the host computer (3) through the stress transmission line (2). Secondly, the ultimate stress value [σ u and the ultimate shear strength [τ f of the fractured rock mass are set in the host computer (3); Finally, according to the real-time change of the stress of the monitored rock mass, the effective stress σ' of the water-containing fissures is analyzed by the water-rock coupling method n and the effective shear stress τ' n , and the degree of fissure cracking is judged by combining fracture mechanics; The ultimate stress value [σ u and the ultimate shear strength [τ f of the fractured rock mass are obtained by taking in-situ samples of the tunnel (6) and using a uniaxial testing instrument to test the uniaxial compressive strength of the rock mass, so as to obtain the compressive strength of the fractured rock mass, which is used as the ultimate stress value [σ u and the ultimate shear strength [τ f , and then the obtained values are input into the host (3) through a mechanical keypad panel.

7. The monitoring and early warning method for tunnel water inrush under the action of water-rock coupling according to claim 6, characterized in that: The prediction of the water inrush risk of rock mass based on machine learning at least includes the following steps: First, collect the geological disaster classification data of the water inrush caused by the tunnel. Collect the changes of the tunnel rock mass stress σ(t) and seepage pressure P(t) in the past half month through the stress sensor (1), collect the crown settlement D(t) at the top of the tunnel in the past half month, and collect the historical maximum rainfall in the past 10 years according to the information of the meteorological bureau for data collection, and further transmit the collected data to the data preprocessing module; Secondly, sort, screen and optimize the data. In the data preprocessing module, use the Z-score method to detect the outliers of the rock mass stress σ(t), seepage pressure P(t) and crown settlement D(t). The expression is: where: Z i is the Z-score of the i-th data point; x i is the value of the i-th data point; μ is the mean of the data set; e is the standard deviation of the data set; Normalize the data at the outlet, and the expression is: In the formula: x is the original data collected; min(x) is the minimum value in the data set; max(x) is the maximum value in the data set; For rainfall feature extraction, extract the statistical features of the historical maximum rainfall R max First, calculate the ratio of the rainfall R(t) in the current period to the historical maximum value Then introduce the attenuation factor λ to weight the influence of historical rainfall: Among them, R(t - k) is the past rainfall with a time step of k; Allocate weights for multiple factors, and construct a data matrix for the index rock mass stress, osmotic pressure, settlement, rainfall ratio, and cumulative rainfall X = [x ij n×m ​ where x ij represents the value of the j-th indicator at the i-th time point, n is the number of time steps, and m is the number of indicators, indicating how many different indicators there are; Calculate the information entropy of the jth index: Among them, is a normalization constant to ensure that the value of information entropy is between [0, 1]; is the normalized probability, representing the relative proportion of the j-th indicator at the i-th time point; P ij is the standardized value of the j-th indicator at the i-th time point; is the sum of the j-th indicator over all time steps; E j is the information entropy, which measures the uncertainty of the j-th indicator; Then calculate the index weight W j , reflecting its contribution degree: In the formula: m is the number of indexes; Immediately build a risk probability model, define the water inrush event label Y ∈ {0, 1}, and the input feature vector is X = [σ norm , P norm , D norm , R ratio , R cum , where σ norm is the normalized value of rock mass stress; P norm is the normalized value of osmotic pressure; D norm is the normalized value of settlement; R ratio is the rainfall ratio; R cum is the cumulative rainfall; The logistic regression model is: Where: Y is the predicted target variable, indicating whether a water inrush disaster occurs. Y = 1 indicates occurrence, and Y = 0 indicates non-occurrence; β0 is the intercept term; βj is the feature coefficient; is a linear combination in the logistic regression model and is transformed into a probability through an exponential function; Finally, construct a risk probability model to predict the water inrush risk of rock mass, and use the maximum likelihood method to construct a log-likelihood function and optimize the solution of parameters: Among them, P i = P(Y = 1|x i ) is the predicted probability of the 0th sample, that is, the water inrush probability calculated by the logistic regression model; y i is the actual label of the ith sample; the parameter estimate value that maximizes lnL(β) is iteratively solved by the Newton-Raphson method When the risk model is constructed, set the early warning dynamic threshold accordingly. First, calibrate the threshold. Based on the historical water inrush event data, determine the optimal threshold θ through the ROC curve, and it is necessary to satisfy: θ = argmax(Youden index) = argmax(sensitivity + specificity - 1) Set the warning dynamic threshold and calculate the risk probability P in real time new = P(Y = 1|x new ), and trigger the warning condition:

Images