A multi-source data risk identification method and system for a blasting disturbance scene

Abstract

The application belongs to the technical field of tunnel engineering geological disaster prevention and control, and discloses a multi-source data risk identification method and system for a blasting disturbance scene; the method comprises the following steps: obtaining structure state data, disturbance source parameters and blasting vibration response data of a target measuring point; determining a filling medium activation sensitive type and matching a critical disturbance threshold value according to the structure state data; screening effective blasting points and predicting disturbance intensity at the target measuring point; calculating an activation risk coefficient based on the predicted disturbance intensity and the critical disturbance threshold value; when the activation risk coefficient exceeds a risk threshold value, combining microseismic activity characteristics to generate a precursor characteristic index, and outputting a sudden gushing water risk identification result; the application establishes a quantitative correlation between geological structure characteristics and blasting disturbance risks, and realizes the pre-judgment and hierarchical control of sudden gushing water risks.

Description

Technical Field

[0001] This application relates to the field of multi-source monitoring data processing technology, and more specifically, to a method and system for multi-source data risk identification in blasting disturbance scenarios. Background Technology

[0002] In disturbed engineering scenarios, risk status identification typically requires the simultaneous use of target area structural status data, disturbance source parameters, and disturbance response data. Target area structural status data characterizes the sensitivity of the target monitoring area under disturbance, disturbance source parameters characterize the location, intensity, and timing of the external disturbance, and disturbance response data characterizes the actual response status of the target monitoring area after the disturbance.

[0003] Existing methods often use the aforementioned data separately for individual judgments. For example, they might determine the static risk level based on the structural state of the target area, estimate the disturbance intensity based on disturbance source parameters, or perform post-event anomaly identification based on disturbance response signals. Because there is a lack of a unified mapping relationship between data from different sources, structural state data is difficult to convert into disturbance thresholds, disturbance source parameters are difficult to correlate with the load-bearing or instability sensitivity of the target area, and anomalous features in the response signals are difficult to differentiate based on the sensitivity type of the target area.

[0004] Taking a blasting disturbance scenario as an example, disturbance source parameters can include the blasting point location, the amount of explosive charge in a single stage, and the initiation delay time. Target area structural state data can include the type of filling medium, the thickness of the insulating layer, and structural integrity data. Disturbance response data can include vibration response or microseismic response data collected after the blast. If judgment is based solely on a single data source, problems such as inconsistencies between static risk and dynamic response, the inability of a uniform threshold to adapt to different structural states, and the inability to accurately attribute abnormal responses can easily arise.

[0005] Patent application CN114723233A discloses a method for predicting tunnel water inrush and mudslide, comprising: step A, acquiring geological data at a predetermined tunnel location; step B, calculating energy loss during seepage based on the acquired geological data; and step C, predicting the rate of water inrush in the tunnel based on the acquired geological data and the energy loss during seepage, which can provide a certain risk reference for tunnel construction. However, the above method still has the following technical problems in practical applications:

[0006] In karst tunnel construction, blasting vibration is one of the main dynamic external disturbances. While existing technologies can identify the geological structure characteristics of concealed karst water bodies and predict peak vibration velocities based on blasting design parameters, a quantitative correlation between the two is lacking. Specifically, geological survey results cannot be translated into blasting vibration control parameters; blasting vibration control standards only set allowable vibration velocities based on the surrounding rock grade, without considering the different responses of different geological structures to vibration. For example, karst caves filled with sandy soil are prone to liquefaction and instability under vibration, while the activation probability of dry, unfilled karst caves is extremely low, yet both are controlled using the same standard in existing technologies. This results in a separation between geological survey information and blasting disturbance control, failing to form a closed-loop predictive chain.

[0007] Therefore, the main problem with existing technology is that it is impossible to predict the risk of hidden karst water bodies being activated by blasting vibrations before blasting operations, and blasting operations become high-risk moments for sudden water inrush without corresponding early warning measures.

[0008] In view of this, this application proposes a multi-source data risk identification method and system for blasting disturbance scenarios to solve the above problems. Summary of the Invention

[0009] To address the problem that existing technologies struggle to uniformly transform target area structural state data, disturbance source parameters, and disturbance response data into calculable risk identification results, leading to reliance on a single threshold and insufficient utilization of dynamic response features for risk state identification under disturbance effects, this application aims to provide a multi-source data risk identification method and system for blasting disturbance scenarios. This system establishes a complete closed-loop control chain from geological structure identification to dynamic early warning, achieving dual protection of static prediction and dynamic verification of sudden water inrush risk, and significantly reducing the risk of sudden water inrush accidents during karst tunnel blasting construction.

[0010] This application provides the following technical solution: a method for multi-source data risk identification in blasting disturbance scenarios, comprising:

[0011] Acquire multi-source input data corresponding to the target measuring point under the blasting disturbance scenario of karst tunnel. The multi-source input data includes the structural state data of the target measuring point, the disturbance source parameters of the blasting operation to be carried out, and the blasting vibration response data.

[0012] Based on the structural state data analysis, the activation sensitivity type corresponding to the filling medium at the target measurement point is obtained, and the corresponding critical disturbance threshold is matched from the preset threshold mapping table according to the activation sensitivity type.

[0013] Based on the disturbance source parameters, the blasting points in the blasting operation to be carried out are screened to obtain the effective blasting points that meet the disturbance influence conditions of the target measuring point. Based on the disturbance source parameters corresponding to the effective blasting points, the predicted disturbance intensity of the blasting operation to be carried out at the target measuring point is predicted.

[0014] By comparing and analyzing the predicted disturbance intensity and the critical disturbance threshold, the activation risk coefficient of the hidden karst water body at the target measuring point is obtained.

[0015] If the activation risk coefficient is greater than the preset risk threshold, microseismic activity characteristics are generated based on the blasting vibration response data, and precursor characteristic indicators are generated based on the activation sensitivity type and microseismic activity characteristics.

[0016] The risk identification results of sudden water inrush corresponding to the target measuring point are generated based on the activation risk coefficient and precursor characteristic indicators.

[0017] A multi-source data risk identification system for blasting disturbance scenarios, comprising implementing the aforementioned multi-source data risk identification method for blasting disturbance scenarios, including:

[0018] Data acquisition module: acquires multi-source input data corresponding to the target measuring point under the blasting disturbance scenario of karst tunnel. The multi-source input data includes the structural state data of the target measuring point, the disturbance source parameters of the blasting operation to be carried out, and the blasting vibration response data.

[0019] Sensitive Classification Module: Based on structural state data analysis, the activation sensitivity type corresponding to the filling medium at the target measurement point is obtained, and the corresponding critical disturbance threshold is matched from the preset threshold mapping table according to the activation sensitivity type;

[0020] Disturbance prediction module: Based on the disturbance source parameters, the blasting points in the blasting operation to be carried out are screened to obtain the effective blasting points that meet the disturbance influence conditions of the target measuring point, and the predicted disturbance intensity of the blasting operation to be carried out at the target measuring point is predicted based on the disturbance source parameters corresponding to the effective blasting points.

[0021] Comparison and Judgment Module: Based on the predicted disturbance intensity and the critical disturbance threshold, a comparative analysis is performed to obtain the activation risk coefficient of the hidden karst water body at the target measuring point;

[0022] Precursor analysis module: If the activation risk coefficient is greater than the preset risk threshold, microseismic activity characteristics are generated based on the blasting vibration response data, and precursor characteristic indicators are generated based on the activation sensitivity type and microseismic activity characteristics.

[0023] Analysis and early warning module: Generates the risk identification results of sudden water inrush corresponding to the target measuring point based on the activation risk coefficient and precursor characteristic indicators.

[0024] This application presents a multi-source data risk identification method and system for blasting disturbance scenarios, highlighting its technical effects and advantages:

[0025] This application does not rely solely on the surrounding rock grade or conventional geological risk level to predict sudden water inrush. Instead, it obtains the activation sensitivity type based on the type of filling medium and the thickness of the impermeable roof at the target measuring point. This allows geological information to be transformed into parameters that reflect the instability sensitivity under blasting vibration, providing a basis for subsequent differentiated critical vibration threshold matching.

[0026] This application matches critical vibration thresholds based on activation sensitivity types, so that hidden karst water bodies at different target measuring points no longer use a uniform blasting control standard, but instead match different critical vibration thresholds based on differences in the thickness of sandy soil filling, cohesive soil filling, no filling, or impermeable roof, thereby improving the pertinence of blasting disturbance risk assessment.

[0027] This application explicitly defines blasting points for blasting operations that are less than or equal to the preset influence distance from the target measuring point as valid blasting points, thus avoiding the use of blasting design parameters of arbitrary blasting points for risk calculation of target measuring points. At the same time, it calculates the predicted peak velocity for multiple target measuring points separately. For cases where the same target measuring point is affected by multiple blasting points, it determines the predicted peak velocity at the corresponding target measuring point according to the maximum value of the segmented predicted peak velocity or the conservative composite result, making the disturbance propagation relationship between the blasting point and the target measuring point clearer.

[0028] This application obtains the activation risk coefficient by predicting the ratio of peak velocity to critical vibration threshold, enabling the geological structure characteristics and blasting disturbance intensity to be compared under the same dimension. This allows for the determination of whether there is a risk of activation of the hidden karst water body at the target measuring point by blasting vibration before blasting operations.

[0029] When the activation risk coefficient is greater than the preset risk threshold, this application further collects blasting vibration response data during the blasting process and matches the precursor judgment threshold according to the activation sensitivity type. This makes the dynamic microseismic precursor identification not simply rely on a uniform threshold, but adapt to the instability sensitivity of the filling medium at the target measuring point, thereby improving the accuracy of identifying the precursor processes of sudden water inrush, such as liquefaction, shear failure, and crack propagation of the filling medium.

[0030] In summary, this application addresses the technical problem that blasting vibrations during karst tunnel construction may activate hidden karst water bodies, leading to sudden water inrushes. This is achieved through a continuous processing chain encompassing activation sensitivity type determination, differentiated threshold matching, vibration peak prediction, activation risk assessment, and microseismic precursor early warning. Existing technologies cannot predict activation risks before blasting operations and lack dynamic precursor identification methods during blasting. This application establishes a quantitative correlation between geological structural characteristics and blasting disturbance risks, realizing a closed-loop chain from static prediction to dynamic verification, and from risk identification to proactive control. This significantly improves the pertinence, scientific rigor, and engineering operability of sudden water inrush prediction. Attached Figure Description

[0031] Figure 1This is a schematic diagram of a multi-source data risk identification method for blasting disturbance scenarios according to this application;

[0032] Figure 2 This is a schematic diagram of the method for obtaining the predicted peak velocity in this application;

[0033] Figure 3 This is a schematic diagram of the method for generating real-time early warning signals based on precursor characteristic indicators in this application;

[0034] Figure 4 This is a schematic diagram of the structure of a multi-source data risk identification system for blasting disturbance scenarios according to this application. Detailed Implementation

[0035] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0036] Example 1

[0037] Please see Figure 1 As shown, this embodiment provides a multi-source data risk identification method for blasting disturbance scenarios, including:

[0038] Acquire multi-source input data corresponding to the target measuring point under the blasting disturbance scenario of karst tunnel. The multi-source input data includes the structural state data of the target measuring point, the disturbance source parameters of the blasting operation to be carried out, and the blasting vibration response data.

[0039] In specific implementations, methods for obtaining multi-source input data corresponding to the target measurement point include:

[0040] The target measuring point is the location of a hidden karst water body, filling medium, or karst anomaly within a preset distance in front of the tunnel face. The preset distance is determined based on the tunnel's single excavation cycle advance, the effective detection distance of advanced geological forecasts, and the impact range of blasting vibration, so that the target measuring point is located within the impact range of the blasting disturbance in front of the tunnel face.

[0041] The structural state data of the target measuring point is used to characterize the sensitive state of the hidden karst water body or filling medium at the target measuring point under the action of blasting disturbance; the structural state data includes drilling data, ground-penetrating radar characteristics and tunnel seismic wave advance prediction characteristics.

[0042] In specific implementations, methods for obtaining drilling data include:

[0043] Advanced drilling holes are laid out at the tunnel face, with the drilling depth covering the target measuring points within a preset distance in front of the tunnel face. During the drilling process, the length of each drilling pass is recorded, and samples of the filling medium, the thickness of the waterproof roof, and core data are extracted as drilling data.

[0044] In specific implementations, methods for obtaining ground-penetrating radar features include:

[0045] A shielded antenna is used to lay transverse survey lines at the tunnel face to obtain radar reflection profiles. The antenna frequency is set based on a balance between detection depth and resolution. At least three transverse survey lines are laid at the tunnel face, located at the distances from the first survey line at the arch, the center of the tunnel face, and the distance from the second survey line at the bottom plate. The distances between the first and second survey lines are determined based on the tunnel cross-section dimensions. Radar data is collected continuously along each survey line, with the antenna bottom pressed against the tunnel face rock wall and the movement speed kept uniform.

[0046] The collected raw radar data were organized into radar gathers according to the survey line number. Each radar gather includes radar channels collected sequentially along the survey line direction. The radar channel numbered m in the radar gather is denoted as Am(n), where m is the radar channel number, n is the sampling point number, and Am(n) is the reflection amplitude of the corresponding sampling point. The radar gathers were then subjected to zero-point correction, DC filtering, gain adjustment, background removal filtering, bandpass filtering, offset repositioning, and time-depth conversion in sequence to obtain the radar reflection profile.

[0047] The zero-point correction method includes: for radar channel Am(n) numbered m, within a preset search sampling interval after the transmission pulse, calculating the absolute value of the difference in reflection amplitude between adjacent sampling points; determining the sampling point whose absolute value first exceeds the jump determination threshold as the jump sampling point of the corresponding radar channel; determining the sampling time corresponding to the jump sampling point as the time zero point of the corresponding radar channel; shifting the sampling time axis of the corresponding radar channel according to the time zero point to align the time zero points of each radar channel; the jump determination threshold is determined based on the average value and standard deviation of the difference in reflection amplitude between adjacent sampling points within the noise sampling interval before the transmission pulse.

[0048] The method for DC component subtraction includes: for radar channel Am(n) numbered m after zero-point correction, calculating the average reflection amplitude of the corresponding radar channel in the sampling interval participating in subsequent processing, and taking the average reflection amplitude as the DC component of the corresponding radar channel; subtracting the corresponding DC component from the reflection amplitude of each sampling point in the corresponding radar channel to obtain the radar channel after DC component subtraction.

[0049] The amplitude compensation method includes: setting an amplitude compensation window that moves along the sampling time direction for the radar channel after DC component subtraction; calculating the square root of the average square of the reflected amplitude within each amplitude compensation window to obtain the window amplitude reference; determining the compensation coefficient based on the ratio of the target amplitude reference to the window amplitude reference; and multiplying the reflected amplitude of the corresponding sampling point by the compensation coefficient to obtain the amplitude-compensated radar channel. The target amplitude reference is taken as the median of the window amplitude references for all radar channels of the corresponding survey line within the sampling interval near the transmitter; the length of the amplitude compensation window is taken as an integer multiple of the number of sampling points corresponding to the duration of the antenna's main pulse.

[0050] The background channel subtraction method includes: in the radar channel set of the same survey line, for the same sampling point number n, calculating the average reflection amplitude of each radar channel participating in the background statistics at the corresponding sampling point number, to obtain the background channel B(n); subtracting B(n) from the reflection amplitude corresponding to sampling point number n in each radar channel Am(n), to obtain the radar channel after background channel subtraction; the radar channels participating in the background statistics are radar channels within the preset number range on both sides of the radar channel to be processed; when the radar channel to be processed is located at the end of the survey line, resulting in insufficient radar channels on one side, radar channels on the other side that meet the preset number of channels are used to participate in the background statistics; the preset number of channels is determined according to the number of radar channels corresponding to the abnormal length threshold, and is not less than the number of radar channels corresponding to the abnormal length threshold.

[0051] The bandpass filtering method includes: determining the lower cutoff frequency fL and the upper cutoff frequency fU based on the nominal center frequency fc and nominal operating bandwidth BW of the radar antenna, where the lower cutoff frequency fL is the difference between fc and BW / 2, and the upper cutoff frequency fU is the sum of fc and BW / 2; performing a Fourier transform on the radar channel after background channel subtraction to obtain the frequency domain signal; retaining the frequency domain signal components with frequencies not less than fL and not greater than fU, and setting the frequency domain signal components with frequencies less than fL and greater than fU to zero; and then performing an inverse Fourier transform to obtain the bandpass filtered radar channel.

[0052] The offset repositioning process includes: performing offset imaging processing on the bandpass-filtered radar gather based on the electromagnetic wave propagation velocity of the medium corresponding to the target survey line, so that the reflected energy corresponding to the diffraction curve is repositioned to the spatial position of the reflector; the electromagnetic wave propagation velocity is obtained through drilling depth calibration; specifically: selecting the anomalous interfaces that have been drilled and exposed within the survey line range, reading the two-way travel time and drilling depth of the corresponding anomalous interface in the radar gather, and obtaining the electromagnetic wave propagation velocity by calculating the ratio of twice the two-way travel time to the drilling depth; when there are more than two drilled and exposed interfaces, calculating the electromagnetic wave propagation velocity corresponding to each drilled and exposed interface separately, and taking the arithmetic mean as the electromagnetic wave propagation velocity of the medium corresponding to the target survey line.

[0053] The time-depth conversion method includes: converting the two-way travel time of the radar reflection signal into depth based on the electromagnetic wave propagation speed; for the reflection signal acquired based on the sampling time, the corresponding depth of the reflection signal is obtained by calculating the product of the electromagnetic wave propagation speed of the medium corresponding to the target survey line and the sampling time, and dividing by 2; after completing the time-depth conversion, a radar reflection profile is generated with the horizontal position of the survey line as the horizontal axis, the depth as the vertical axis, and the reflection amplitude after time-depth conversion as the pixel value.

[0054] The first anomalous region is identified from the radar reflection profile. The identification method is as follows: Regions where the ratio of the reflected wave amplitude to the average reflected wave amplitude of the background region is greater than or equal to an amplitude multiplier threshold, and the continuous reflection length is not less than anomaly length threshold, are marked as the first anomalous region. Specifically, the reflected wave amplitude is obtained by extracting the reflected wave amplitude values ​​of each data track from the processed radar reflection profile. For each data track in the radar reflection profile, the maximum amplitude value of the reflected wave within the corresponding data time window is taken, i.e., the maximum absolute value from the wave crest to the wave trough, as the reflected wave amplitude of the corresponding track. The background region is selected by determining a normal rock mass area as the background region along the survey line direction. The background region should be selected where the coefficient of variation of the reflected wave amplitude in the radar image is less than the background coefficient of variation threshold, and no continuous strong reflection interface appears. The corresponding segment should be located within a certain distance range on both sides of the potential anomalous region. If the anomalous region is located at the end of the survey line, only one side of the background region may be selected.

[0055] Calculation of the average reflection amplitude of the background area: In the selected background area, extract the main amplitude value of each channel along the survey line direction, calculate the arithmetic mean of all extracted amplitude values, and obtain the average reflection amplitude of the background area.

[0056] Method for setting the amplitude multiplier threshold: The amplitude multiplier threshold is used to distinguish between the reflection of anomalies and the reflection of background rock masses. Karst anomalies, especially water-bearing or filled karst caves, have significant differences in wave impedance compared with intact surrounding rock. The amplitude multiplier threshold can be adjusted according to the reflection characteristics of known geological calibration points on site. The calibration method is as follows: During tunnel excavation, select the location of typical karst anomalies that have been exposed, compare the radar reflection profile obtained by advance prediction with the actual geological conditions at the corresponding location, and statistically analyze the ratio distribution of the anomaly reflection amplitude to the average reflection amplitude of the background area. Take the lower limit of the ratio distribution as the amplitude multiplier threshold.

[0057] The background area is selected as follows: In the radar reflection profile, the initial radar trace number of the candidate anomaly area along the survey line direction is denoted as... The radar channel number at the end will be recorded as The number of radar channels contained in the candidate anomaly region is denoted as , ; number of intervals Set as Round up to the nearest integer value to determine the background channel count. Set as ;Will to The corresponding radar channel is used as the candidate background area on the left. to The corresponding radar channel is used as the right candidate background area; when the left or right candidate background area exceeds the boundary of the survey line, the candidate background area that can be obtained within the boundary of the survey line is used as the background area; when candidate background areas can be obtained on both the left and right sides, the candidate background areas on both sides are used together as the background area.

[0058] The method for calculating the average reflection amplitude of the background region is as follows: In the background region, radar reflection data within the same depth range as the candidate anomaly region are selected; for each radar channel in the background region, the maximum absolute value of the reflection amplitude of the corresponding radar channel within the depth range is calculated as the main reflection amplitude of the corresponding radar channel; the arithmetic mean of the main reflection amplitudes of all radar channels in the background region is calculated to obtain the average reflection amplitude of the background region.

[0059] The method for calculating the reflection amplitude ratio of the candidate anomaly region is as follows: In the candidate anomaly region, radar reflection data within the same depth range as the background region are selected; for each radar channel in the candidate anomaly region, the maximum absolute value of the reflection amplitude of the corresponding radar channel within the depth range is calculated; the maximum value of the main reflection amplitude of all radar channels in the candidate anomaly region is divided by the average reflection amplitude of the background region to obtain the reflection amplitude ratio.

[0060] The reflection amplitude threshold is the product of the amplitude multiplier threshold and the average reflection amplitude of the background area; the method for setting the anomaly length threshold is: determined according to the tunnel cross-section size, detection accuracy requirements and anomaly size characteristics. The calibration method is: during tunnel excavation, select the location of typical exposed karst anomalies, count the continuous reflection length of the anomalies in the radar reflection profile, and take the lower limit of the corresponding length as the anomaly length threshold.

[0061] The ratio of the maximum reflection amplitude in the first anomaly region to the average reflection amplitude in the background region is calculated to obtain the reflection amplitude ratio. The average reflection amplitude in the background region is the arithmetic mean of the principal amplitude values ​​in the background region. The reflection amplitude ratio is used as the output ground-penetrating radar feature for subsequent joint determination of activation sensitive types.

[0062] In specific implementations, methods for obtaining tunnel seismic wave advance prediction characteristics include:

[0063] Receiver holes and boreholes are installed on both sides of the tunnel sidewalls. Three-component geophones are installed in the receiver holes, each capable of simultaneously recording three orthogonal components of the seismic wave propagation process: the horizontal transverse component, the horizontal longitudinal component, and the vertical component. The longitudinal wave energy is mainly concentrated in the component perpendicular to the propagation direction, while the transverse wave energy is mainly concentrated in the component perpendicular to the propagation direction. The geometric parameters of the observation system, including receiver hole depth, borehole depth, receiver hole spacing, borehole spacing, minimum offset, number of boreholes, and number of receiver holes, are determined based on the detection distance, surrounding rock conditions, and resolution requirements to ensure that the seismic wave signal covers the entire detection range without signal overlap. For example, the receiver hole layout method is as follows: the number of receiver holes on each sidewall is determined according to the detection coverage range to ensure that seismic wave signals from different directions in front of the tunnel face can be received; the receiver hole depth is determined based on the required depth to ensure the detection... The principles for determining the parameters are: good coupling between the detector and the surrounding rock, and avoiding the loosening zone of the excavation; the receiver hole diameter is determined based on the outer diameter of the geophone to ensure that the geophone can be smoothly inserted and coupled with the hole wall; the receiver hole spacing is determined based on the principle of ensuring that the spatial sampling interval meets the wave field reconstruction requirements; the borehole layout method is as follows: boreholes are laid on each side wall behind the receiver holes, and the number of boreholes is determined according to the detection distance and excitation energy requirements to ensure that the seismic wave signal can cover the entire detection range; the borehole depth is determined based on the principle of ensuring effective propagation of excitation energy without damaging the receiver holes; the borehole diameter is determined based on the diameter of the explosive charge to ensure that the explosive can be smoothly inserted; the borehole spacing is determined based on the principle of ensuring that the signals of adjacent boreholes do not overlap and cover the entire detection section; boreholes are detonated one by one to obtain seismic wave propagation data; the specific values ​​of the above parameters can be adjusted according to the tunnel surrounding rock conditions, the model of the tunnel seismic wave advance prediction system, and the results of field tests.

[0064] Detonation is performed hole by hole to generate seismic waves. As these waves propagate through the rock mass, they are reflected when they encounter interfaces with different wave impedances, such as the boundaries of karst anomalies. These reflected waves are received by geophones, which record the seismic wave waveform data of the horizontal transverse component, the horizontal longitudinal component, and the vertical component acquired by each receiver channel. Based on the strength of the waveform amplitude in each component and the polarization characteristics of the particle vibration, P-waves and S-waves can be distinguished: P-waves have particle vibration in the same direction as the propagation direction, and the amplitude is strongest in the component with the wave propagation direction; S-waves have particle vibration in the direction perpendicular to the propagation direction, and the amplitude is strongest in the component perpendicular to the propagation direction.

[0065] The acquired seismic wave data is imported into a pre-set data processing software for preprocessing and then separated from the P-wave and S-wave data to obtain P-wave and S-wave data volumes.

[0066] First-break picking is performed on the shear wave data volume, that is, identifying the time point when the shear wave arrives at each receiver. The first-break picking methods include but are not limited to the ratio of the mean of long-time window to short-time window method or manual interactive picking method. After picking, visual inspection is performed and manual correction is performed on first-break points that deviate significantly, ensuring that the picking error is controlled within one sampling interval.

[0067] According to the straight-line spatial distance between the blast hole and the receiver, as well as the first-break time of the shear wave, the average shear wave velocity of the rock mass on the corresponding ray path is calculated, which is the ratio of the straight-line spatial distance to the first-break time of the shear wave. The calculation is performed for every combination of each blast hole and each receiver to obtain shear wave velocity values at discrete points covering the detection area. The average P-wave velocity is calculated by the same method and will not be repeated in detail.

[0068] Statistics or inversion are performed on the P-wave velocity and Poisson's ratio calculated from different combinations of blast holes and receivers in the exploration section to obtain the spatial distribution of P-wave velocity and the spatial distribution of Poisson's ratio in the detection area; specifically, the detection area is discretized into grid cells, with the P-wave velocity values calculated from combinations of blast holes and receivers as input, an inversion algorithm is used to reconstruct the P-wave velocity value in each grid cell to obtain the spatial distribution of P-wave velocity. Specifically, theoretical travel time is calculated through ray tracing, travel time residuals are obtained by comparing the difference between the theoretical travel time and the measured first-break travel time, the velocity model is iteratively updated with the goal of minimizing the travel time residual until the convergence condition is satisfied, and the spatial distribution image of P-wave velocity is output. The setting method of the convergence condition is: it can be adjusted according to the tunnel surrounding rock conditions, the model of the tunnel seismic wave advanced prediction system and the field test results; the specific calculation process of the above inversion method can refer to the operating specifications of the data processing software matched with the tunnel seismic wave advanced prediction system or the common knowledge in the field of seismic tomography.

[0069] Scanning in the detection area, a connected region where the absolute value of the velocity anomaly rate of P-wave velocity relative to the velocity of the background region is not less than the velocity anomaly rate threshold and the continuous spatial length is not less than the anomaly length threshold is marked as a second abnormal region; the setting method of the velocity anomaly rate threshold is: it is determined according to the velocity decrease characteristic of karst fillers relative to intact surrounding rock, and the specific value can be adjusted according to the variation coefficient of the background velocity of the surrounding rock and the signal-to-noise ratio; during tunnel excavation, positions of known typical karst abnormal bodies can be selected, the distribution range of the velocity anomaly rate of the abnormal bodies is counted, and the lower limit value of the distribution range of the velocity anomaly rate of the abnormal bodies is taken as the velocity anomaly rate threshold; the setting method of the anomaly length threshold is: it is determined according to the spatial extension size of the karst abnormal body, known typical karst abnormal bodies can be selected during excavation, their continuous reflection length presented in tunnel seismic wave advanced prediction detection is counted, and the lower limit value of the continuous reflection length is taken as the anomaly length threshold.

[0070] The difference between the P-wave velocity in the anomalous area and the P-wave velocity in the background area is calculated as a percentage of the P-wave velocity in the background area to obtain the P-wave velocity anomaly rate. The background area is the tunnel seismic wave advance prediction data within the background distance range in front of the anomalous area. The method for selecting the background area is as follows: with the candidate anomalous area as a reference, an area within a certain distance range in front of the candidate anomalous area is selected as the background area along the tunnel excavation direction. The background distance range is determined to ensure that the background wave velocity statistics are stable and not affected by the anomalous body.

[0071] Poisson's ratio is calculated based on the P-wave velocity and S-wave velocity. Specifically, Poisson's ratio... ,in, For longitudinal wave velocity; The shear wave velocity is given; for the second anomaly region, the average value of the Poisson's ratio calculated from all target measurement points within the second anomaly region is taken as the output.

[0072] The P-wave velocity anomaly rate and Poisson's ratio are used as output tunnel seismic wave advance prediction features for subsequent joint determination of activation sensitive types.

[0073] The disturbance source parameters for the blasting operation to be carried out are used to characterize the conditions under which the blasting disturbance source acts on the target measuring point; the disturbance source parameters include the location of the blasting point, the amount of explosive charge in a single stage, and the detonation delay time.

[0074] Blasting vibration response data is used to characterize the dynamic response state of the target measuring point under blasting disturbance. Blasting vibration response data includes microseismic data collected before and after blasting, and the microseismic data includes the occurrence time, energy density, and spectral distribution of microseismic events.

[0075] In specific implementations, methods for acquiring microseismic data include:

[0076] From the collected microseismic data, individual microseismic events are identified. The method for determining a microseismic event is as follows: when the signal amplitude exceeds the event trigger threshold, and the duration of continuous exceedance of the threshold is not less than the minimum duration of the event, it is determined to be a microseismic event. The method for setting the event trigger threshold is: dynamically set according to the statistical characteristics of background noise. The method for setting the minimum duration of the event is: collect background noise data when there are no effective microseismic events at the tunnel construction site, statistically analyze the duration distribution of instantaneous impulse noise in the background noise, and take the upper limit of the impulse noise duration as the minimum duration of the event.

[0077] For each microseismic event, the moment when the signal amplitude first exceeds the event trigger threshold is taken as the occurrence time; a fast Fourier transform is performed on the signal within the corresponding time window of the microseismic event to obtain the power spectrum; the power spectrum is integrated to obtain the total energy, and the total energy is divided by the time window length to obtain the energy density; the frequency corresponding to the maximum power spectral density is taken as the dominant frequency of the microseismic event.

[0078] The setting method of the time window is as follows: taking the occurrence time as the starting point, the event duration is taken backward; the event duration is the time interval from when the signal amplitude first exceeds the trigger threshold to when it last drops below the trigger threshold, and an additional event extension duration is added. The setting method of the event extension duration is as follows: collect typical microseismic event waveforms at the tunnel construction site, analyze the time required for the microseismic event to attenuate from the amplitude below the trigger threshold to fully recover to the background noise level, take the statistical average as the event extension duration, or set it according to the default reference value of the data processing software supporting the microseismic monitoring system.

[0079] Calculate background parameters and characteristic parameters respectively according to microseismic data collected within a first preset time before blasting and within a second preset time after blasting; the background parameters include a background value of microseismic event frequency and a background value of microseismic energy density; the background value of microseismic event frequency is the number of microseismic events recorded within the first preset time before blasting; the background value of microseismic energy density is the average energy density of all microseismic events within the first preset time before blasting; the characteristic parameters include microseismic event frequency, microseismic energy density and proportion of high-frequency microseismic signals; wherein the microseismic event frequency is the number of microseismic events recorded within the second preset time after blasting; the microseismic energy density is the average energy density of all microseismic events within the second preset time after blasting; the proportion of high-frequency microseismic signals is the percentage of the number of microseismic events with main frequency higher than a frequency threshold in the total number of microseismic events within the second preset time after blasting. The setting method of the frequency threshold is as follows: collect long-period background noise data during periods without blasting operation and significant construction operation interference at the tunnel construction site, perform fast Fourier transform on the background noise to obtain spectral distribution, calculate the 1 / 3 octave spectrum of the power spectral density of the background noise, set the frequency threshold near the starting frequency where the power spectral density of the background noise drops to a flat section, and can be adjusted according to actual conditions.

[0080] Calculate precursory characteristic indicators according to the background parameters and the characteristic parameters, wherein the precursory characteristic indicators include microseismic frequency anomaly rate, energy density anomaly rate and proportion of high-frequency microseismic signals. Specifically, calculate the ratio of the difference between the microseismic event frequency and the background value of microseismic event frequency to the background value of microseismic event frequency, to obtain the microseismic frequency anomaly rate; calculate the ratio of the difference between the microseismic energy density and the background value of microseismic energy density to the background value of microseismic energy density, to obtain the energy density anomaly rate.

[0081] Obtain the activation sensitivity type corresponding to the filling medium at the target measuring point according to the structural state data, and match the corresponding critical disturbance threshold from a preset threshold mapping table according to the activation sensitivity type;

[0082] In specific embodiments, the method for obtaining the activation sensitivity type corresponding to the filling medium at the target measuring point comprises:

[0083] The type of filling medium is determined based on the results of particle size distribution test and liquid limit and plastic limit test. Specifically, the filling medium sample is dried and sieved to obtain the total dry mass of the filling medium sample. The sieve aperture for sand division is the sieve aperture used in the particle size distribution test to separate the sand group and the silt group.

[0084] The method for setting the sand subdivision boundary sieve aperture is as follows: Before conducting the particle analysis test of the filling medium, determine the soil particle group division basis used in this test, and read the particle size boundary value between the sand particle group and the silt particle group from the particle group division basis. If there is a standard sieve in the sieve group used in this particle analysis test with a sieve aperture diameter equal to the particle size boundary value, then the standard sieve is determined as the sand subdivision boundary sieve aperture.

[0085] If no standard sieve with an aperture diameter equal to the particle size boundary value exists in the sieve group used in this particle size analysis experiment, the particle size boundary value will be used as the equivalent sand size boundary sieve aperture diameter, and the mass passing through the particle size boundary value will be calculated based on the particle size distribution curve. Specifically, the sieve aperture diameters will be arranged in descending order to determine the sieve aperture diameters that are larger than the particle size boundary value and adjacent to it. And sieve pore sizes smaller than the particle size boundary value and adjacent to the particle size boundary value. Read them separately Corresponding cumulative pass rate and Corresponding cumulative pass rate On a particle size distribution curve where the horizontal axis represents the logarithm of the sieve aperture and the vertical axis represents the cumulative passing mass percentage, linear interpolation is performed to calculate the cumulative passing mass percentage corresponding to the particle size boundary value. In one embodiment, the cumulative passing mass percentage corresponding to the particle size boundary value is... The calculation formula is: ;in, This is the particle size cutoff value.

[0086] The dry mass of the fine-grained group with a particle size smaller than the particle size boundary value is obtained by multiplying the total dry mass of the filling medium sample with the cumulative mass corresponding to the particle size boundary value by the mass percentage. The dry mass of the sand group with a particle size greater than or equal to the particle size boundary value is the difference between the total dry mass of the filling medium sample and the dry mass of the fine-grained group.

[0087] If the particle analysis test report has already given the dry mass results of the sand group, silt group and clay group, then the dry mass of the sand group and the dry mass of the fine group will no longer be calculated by interpolation of the sieve size. Instead, the dry mass of the sand group will be read directly, and the sum of the dry mass of the silt group and the dry mass of the clay group will be used as the dry mass of the fine group.

[0088] Liquid limit and plastic limit tests were performed on the fine particles passing through the sand fine boundary sieve to obtain the liquid limit and plastic limit, and the difference between the liquid limit and plastic limit was calculated to obtain the plasticity index.

[0089] If the mass ratio is greater than or equal to the plasticity threshold and the plasticity index is less than 7, the filling medium type is determined to be sandy soil; if the mass ratio is less than the plasticity threshold and the plasticity index is greater than or equal to 7, the filling medium type is determined to be cohesive soil; if the drilling record shows that no effective filling medium was obtained inside the anomaly, and the geological radar characteristics and tunnel seismic wave advance prediction characteristics do not meet the filling anomaly determination conditions, the filling medium type is determined to be unfilled.

[0090] The method for setting the plasticity judgment threshold is as follows: Perform at least two parallel liquid limit tests and at least two parallel plastic limit tests on the same filling medium sample, obtaining the maximum and minimum liquid limit values ​​from the parallel liquid limit test results, and the maximum and minimum plastic limit values ​​from the parallel plastic limit test results. The sum of the ranges from the parallel liquid limit tests and the parallel plastic limit tests is used as the plasticity judgment threshold. The range from the parallel liquid limit tests is the difference between the maximum and minimum liquid limit values, and the range from the parallel plastic limit tests is the difference between the maximum and minimum plastic limit values. When the plasticity index is greater than the plasticity judgment threshold, the fine-particle portion is determined to have identifiable plasticity; when the plasticity index is less than or equal to the plasticity judgment threshold, the fine-particle portion is determined to lack identifiable plasticity.

[0091] The method for determining the type of filling medium is as follows: When the dry mass of coarse particles is greater than that of fine particles, and the fine particles do not exhibit identifiable plasticity, the filling medium type is determined to be sandy soil; when the dry mass of fine particles is greater than or equal to that of coarse particles, and the fine particles exhibit identifiable plasticity, the filling medium type is determined to be cohesive soil; when the dry mass of coarse particles is greater than that of fine particles, but the fine particles exhibit identifiable plasticity, the influence of fine particle plasticity is calculated, and the corresponding influence of fine particle plasticity is compared with the plasticity threshold. If the influence of fine particle plasticity is greater than the plasticity threshold, the filling medium type is determined to be cohesive soil; otherwise, it is determined to be sandy soil; when the dry mass of fine particles is greater than or equal to that of coarse particles, but the fine particles do not exhibit identifiable plasticity, the filling medium type is determined to be sandy soil; in one embodiment, the influence of fine particle plasticity... The calculation formula is: ;in, This refers to the dry weight of coarse particles; Dry weight of fine particles; It is the plasticity index.

[0092] If the drilling record shows that no effective filling medium sample was obtained inside the anomaly, it is not directly determined to be without filling; further judgment is made on whether there is an abnormal filling response based on the characteristics of ground-penetrating radar and the characteristics of tunnel seismic wave advance prediction.

[0093] When the characteristics of the ground-penetrating radar do not meet the criteria for determining radar filling anomalies, and the characteristics of tunnel seismic wave advance prediction do not meet the criteria for determining low-velocity filling anomalies, the filling medium type is determined to be unfilled.

[0094] When the characteristics of the ground-penetrating radar meet the criteria for radar filling anomaly determination, or the characteristics of tunnel seismic wave advance prediction meet the criteria for low-velocity filling anomaly determination, the type of filling medium cannot be determined by sample detection, and the process enters the joint determination process of core quality index value, ground-penetrating radar characteristics and tunnel seismic wave advance prediction characteristics.

[0095] The criteria for determining radar filling anomalies are as follows: the reflection amplitude ratio of the candidate anomaly region is greater than the upper limit of radar background fluctuation; the reflection amplitude ratio of the candidate anomaly region is the ratio of the main reflection amplitude of the candidate anomaly region to the average reflection amplitude of the background region; the method for setting the upper limit of radar background fluctuation is as follows: within the background region, according to the same number of radar channels and the same depth range as the candidate anomaly region, background sub-regions are sequentially intercepted along the survey line direction; the ratio of the main reflection amplitude of each background sub-region to the average reflection amplitude of the background region is calculated respectively; the maximum value among the corresponding ratios of each background sub-region is taken as the upper limit of radar background fluctuation; if the reflection amplitude ratio of the candidate anomaly region is greater than the upper limit of radar background fluctuation, the geological radar feature is determined to meet the criteria for determining radar filling anomalies; otherwise, the geological radar feature is determined not to meet the criteria for determining radar filling anomalies.

[0096] The criteria for determining low-speed filling anomalies are: the longitudinal wave velocity in the candidate anomaly region is less than the longitudinal wave velocity in the background region, and the anomaly rate of the longitudinal wave velocity in the candidate anomaly region is greater than the upper limit of the background wave velocity fluctuation.

[0097] The P-wave velocity anomaly rate of the candidate anomaly region is obtained by dividing the difference between the P-wave velocity in the background region and the P-wave velocity in the candidate anomaly region by the P-wave velocity in the background region; the P-wave velocity in the background region is taken as the arithmetic mean of the P-wave velocities in the background region, and the P-wave velocity in the candidate anomaly region is taken as the arithmetic mean of the P-wave velocities in the candidate anomaly region.

[0098] The method for setting the upper limit of background wave velocity fluctuation is as follows: Within the background region, background sub-regions are sequentially extracted according to the same number of grid cells and the same spatial shape as the candidate anomaly region; the wave velocity deviation rate of each background sub-region relative to the P-wave velocity of the background region is calculated, and the wave velocity deviation rate is the ratio of the difference between the P-wave velocity of the background region and the average P-wave velocity of the background sub-region to the P-wave velocity of the background region; the maximum value of the wave velocity deviation rate corresponding to each background sub-region is taken as the upper limit of background wave velocity fluctuation.

[0099] If the P-wave velocity in the candidate anomaly region is less than that in the background region, and the anomaly rate of the P-wave velocity in the candidate anomaly region is greater than the upper limit of the background velocity fluctuation, then the tunnel seismic wave advance prediction characteristics are determined to meet the low-velocity filling anomaly judgment conditions; otherwise, the tunnel seismic wave advance prediction characteristics are determined not to meet the low-velocity filling anomaly judgment conditions.

[0100] The thickness of the water-tight roof is read from the top of the abnormal body in the borehole record, such as the drilling speed change point or the starting point of circulating fluid leakage, and the vertical distance to the working face. When there are multiple vertical distances, each vertical distance is read separately, and the minimum value among the multiple vertical distances is taken as the thickness of the water-tight roof. The core data is the cumulative length of the cores in each run whose length is greater than or equal to the core length threshold, and the corresponding run footage.

[0101] The criteria for determining drilling data include:

[0102] If the filling medium is sandy soil, or the thickness of the waterproof roof is less than the first thickness threshold, then the activation sensitivity type is determined to be Class A, corresponding to the high sensitivity type;

[0103] Otherwise, if the filling medium is cohesive soil, or the thickness of the waterproof roof is between the first thickness threshold and the second thickness threshold, including the first thickness threshold and the second thickness threshold, then the activation sensitivity type is determined to be Class B, corresponding to the medium sensitivity type.

[0104] Otherwise, if the filling medium type is unfilled, or the thickness of the waterproof roof is greater than the second thickness threshold, the activation sensitivity type is determined to be Class C, corresponding to the low sensitivity type; among them, the activation sensitivity of Class A, Class B and Class C decreases in that order; when the determination results corresponding to the filling medium type and the waterproof roof thickness are inconsistent, the activation sensitivity type with the higher risk level is taken as the final activation sensitivity type; the first thickness threshold corresponds to the minimum safe thickness that can effectively block the transmission of blasting vibration energy, and the second thickness threshold corresponds to the thickness limit that significantly reduces the risk of instability of the filling medium. The method for setting the first thickness threshold and the second thickness threshold is: it can be calibrated by numerical simulation or field test based on the thickness of the overburden layer at the top of the tunnel.

[0105] If at least one of the filling medium type and the thickness of the impermeable roof cannot be determined, or if the two determinations are inconsistent, the core quality index value is obtained based on the core data analysis. The core quality index value is specifically the ratio of the cumulative length of cores with a length greater than or equal to the core length threshold in each core run to the corresponding run's advance length. Ground-penetrating radar (GPR) data and tunnel seismic wave advance prediction data are collected, and GPR characteristics and tunnel seismic wave advance prediction characteristics are obtained based on the analysis of the GPR data and tunnel seismic wave advance prediction data. The activation sensitivity type is jointly determined by combining the core quality index value, GPR characteristics, and tunnel seismic wave advance prediction characteristics. The core length threshold is set according to the relevant rock mass classification standards, corresponding to the minimum core length that can reflect the integrity of the rock mass.

[0106] If the core quality index value is less than the first core quality index threshold, and the Poisson's ratio in the tunnel seismic wave advance prediction data is not less than the Poisson's ratio threshold, and the reflection amplitude ratio in the ground-penetrating radar data is not less than the amplitude ratio threshold, then the activation sensitivity type is determined to be Class A.

[0107] If the core quality index value is between the first core quality index threshold and the second core quality index threshold, and the longitudinal wave velocity anomaly rate in the tunnel seismic wave advance prediction data is between the first wave velocity anomaly rate threshold and the second wave velocity anomaly rate threshold, then the activation sensitivity type is determined to be Class B.

[0108] If the core quality index value is greater than the second core quality index threshold and the longitudinal wave velocity anomaly rate is less than the first wave velocity anomaly rate threshold, then the activation sensitivity type is determined to be Class C.

[0109] Otherwise, a default judgment is made based on the core quality index value: a core quality index value less than the first core quality index threshold is classified as Class A; a core quality index value between the first and second core quality index thresholds is classified as Class B; and a core quality index value greater than the second core quality index threshold is classified as Class C. The method for setting each threshold is as follows: the first and second core quality index thresholds are determined based on the core quality index values ​​and activation sensitivity type labels in the calibration sample set. Specifically, each sample in the calibration sample set is sorted in ascending order of core quality index value. When the activation sensitivity type labels of two adjacent samples are inconsistent, the arithmetic mean of the core quality index values ​​of the corresponding two adjacent samples is used as the candidate core quality index boundary point. Two candidate boundary points are selected from the candidate core quality index boundary points and used as the first and second candidate core quality index thresholds, respectively. The calibration sample set is then classified according to the following rules: when the core quality index value of a sample... If the core quality index value is less than the threshold of the first candidate core quality index, the corresponding sample is classified as Class A; if the core quality index value of the sample is greater than or equal to the threshold of the first candidate core quality index and less than the threshold of the second candidate core quality index, the corresponding sample is classified as Class B; if the core quality index value of the sample is greater than or equal to the threshold of the second candidate core quality index, the corresponding sample is classified as Class C. The number of inconsistent samples between the classification result corresponding to each group of candidate boundary points and the activated sensitive type label in the calibrated sample set is calculated one by one. The group of candidate boundary points with the smallest number of inconsistent samples is determined as the first core quality index threshold and the second core quality index threshold. When there are two or more groups of candidate boundary points with the same number of inconsistent samples, the sum of the distances from each group of candidate boundary points to the core quality index values ​​of the nearest samples on both sides is calculated. The group of candidate boundary points with the largest sum of distances is determined as the first core quality index threshold and the second core quality index threshold. If the sum of distances is still the same, the group of candidate boundary points that minimizes the number of missed Class A samples is selected.

[0110] The Poisson ratio threshold is determined based on the Poisson ratio μ and activation sensitivity type labels in the calibration sample set. Specifically, the samples in the calibration sample set are sorted in ascending order of their Poisson ratio μ. When two adjacent samples have different activation sensitivity type labels (one being class A and the other not), the arithmetic mean of the Poisson ratios of the two adjacent samples is used as a candidate Poisson ratio cutoff point. Each candidate Poisson ratio cutoff point is used as a candidate Poisson ratio threshold, and the calibration sample set is classified into two classes according to the following rule: when a sample's Poisson ratio is greater than or equal to the candidate Poisson ratio threshold, the corresponding sample is classified as full. The sample must meet the Poisson's ratio condition for Class A. When the sample Poisson's ratio is less than the candidate Poisson's ratio threshold, the corresponding sample is judged as not meeting the Class A Poisson's ratio condition. The number of inconsistent samples between the judgment result corresponding to each candidate Poisson's ratio threshold and the Class A label in the calibration sample set is calculated, and the candidate Poisson's ratio threshold with the smallest number of inconsistent samples is determined as the Poisson's ratio threshold. When there are two or more candidate Poisson's ratio thresholds with the same number of inconsistent samples, the candidate Poisson's ratio threshold with the smallest number of missed Class A samples is selected first. If the number of missed Class A samples is still the same, the candidate Poisson's ratio threshold with the largest sum of the Poisson's ratio distances to the nearest samples on both sides is selected.

[0111] The amplitude ratio threshold is determined based on the reflection amplitude ratio and activation sensitivity type label in the calibration sample set. Specifically, the samples in the calibration sample set are sorted in ascending order of reflection amplitude ratio. When two adjacent samples have different activation sensitivity type labels (e.g., A and B), the arithmetic mean of the reflection amplitude ratios of the two adjacent samples is used as the candidate amplitude ratio cutoff point. Each candidate amplitude ratio cutoff point is used as a candidate amplitude ratio threshold, and the calibration sample set is binary classified according to the following rule: when the reflection amplitude ratio of a sample is greater than or equal to the candidate amplitude ratio threshold, the corresponding sample is classified as full-range. The system assumes that the radar reflection condition meets Class A criteria. When the sample reflection amplitude ratio is less than the candidate amplitude ratio threshold, the corresponding sample is judged as not meeting the Class A radar reflection condition. The system calculates the number of inconsistent samples between the judgment result corresponding to each candidate amplitude ratio threshold and the Class A label in the calibration sample set, and determines the candidate amplitude ratio threshold with the smallest number of inconsistent samples as the amplitude ratio threshold. When there are two or more candidate amplitude ratio thresholds with the same number of inconsistent samples, the candidate amplitude ratio threshold with the smallest number of missed Class A samples is selected first. If the number of missed Class A samples is still the same, the candidate amplitude ratio threshold with the largest sum of the reflection amplitude ratio distances to the nearest samples on both sides is selected.

[0112] The first and second wave velocity anomaly rate thresholds are determined based on the P-wave velocity anomaly rate and activation sensitivity type label in the calibrated sample set. Specifically, the samples in the calibrated sample set are sorted in ascending order of P-wave velocity anomaly rate. When the activation sensitivity type label of two adjacent samples is inconsistent, the arithmetic mean of the P-wave velocity anomaly rates of the corresponding two adjacent samples is used as the candidate wave velocity anomaly rate cutoff point. Two candidate cutoff points are selected from the candidate wave velocity anomaly rate cutoff points and used as the first and second candidate wave velocity anomaly rate thresholds, respectively. The calibrated sample set is then classified according to the following rules: when the P-wave velocity anomaly rate of a sample is less than the first candidate wave velocity anomaly rate threshold, the corresponding sample is classified as class C; when the P-wave velocity anomaly rate of a sample is greater than or equal to the first candidate wave velocity anomaly rate threshold and less than the second candidate wave velocity anomaly rate threshold, the sample is classified as class C. When the candidate wave velocity anomaly rate threshold is reached, the corresponding sample is classified as Class B. When the sample P-wave velocity anomaly rate is greater than or equal to the second candidate wave velocity anomaly rate threshold, the corresponding sample is classified as Class A. The number of inconsistent samples between the classification result corresponding to each group of candidate boundary points and the activation sensitive type label in the labeled sample set is calculated one by one. The group of candidate boundary points with the smallest number of inconsistent samples is determined as the first wave velocity anomaly rate threshold and the second wave velocity anomaly rate threshold. When there are two or more groups of candidate boundary points with the same number of inconsistent samples, the sum of the distances from each group of candidate boundary points to the P-wave velocity anomaly rates of the nearest samples on both sides is calculated. The group of candidate boundary points with the largest sum of distances is determined as the first wave velocity anomaly rate threshold and the second wave velocity anomaly rate threshold. If the sum of distances is still the same, the group of candidate boundary points that minimizes the number of missed Class A samples is selected.

[0113] Based on the activation sensitivity type, the corresponding critical vibration threshold is matched according to the preset threshold mapping table;

[0114] In specific implementations, the method for matching the corresponding critical vibration threshold includes:

[0115] The corresponding critical vibration threshold is matched from a preset threshold mapping table based on the activation sensitivity type; the preset threshold mapping table includes the critical vibration threshold corresponding to each activation sensitivity type; the method for establishing the preset threshold mapping table includes:

[0116] For type A, which corresponds to a highly sensitive type, the instability mode is liquefaction of sandy soil or sliding along the top surface under vibration. The critical vibration threshold is calculated based on Newton's second law and the limit equilibrium method. The critical vibration threshold is calculated when the sensitivity type is activated as type A. ,in, To activate the cohesive force of the filler when the sensitivity type is A, To activate the rock mass integrity coefficient when the sensitivity type is A, To activate the natural density of the infill material when the sensitivity type is A, it can be obtained by taking a core sample and measuring it using the ring cutter method or the wax sealing method. This is a conventional technique and will not be described in detail here. For longitudinal wave velocity; It is the acceleration due to gravity; The thickness of the impermeable roof when the sensitivity type is A is obtained from the advanced horizontal drilling record; To activate the density of the infill material when the sensitivity type is A, it can be obtained by taking a core sample and measuring it using the ring cutter method or the wax sealing method. This is a conventional technique and will not be described in detail here. To activate the internal friction angle of the infill material when the sensitivity type is A, it was obtained through experimental calibration. Specifically, samples of infill material with sandy soil were taken, and R1 groups of different vertical pressures were set. Each group had at least k parallel tests. Horizontal shear force was applied at a preset shear rate until the sample infill material was sheared. The vertical pressure and the corresponding maximum shear stress for each group were recorded. A straight line corresponding to the maximum shear stress was fitted with respect to the vertical pressure, and the slope of the line was taken as the tangent of the internal friction angle. The internal friction angle was calculated based on the tangent of the internal friction angle. The method for setting R1 is as follows: the vertical pressure range is determined according to the actual stress range that the infill material may bear in the formation, and the number of groups is set according to the vertical pressure range. ,like ,in, This represents the upper limit of the vertical pressure range. This represents the lower limit of the vertical pressure range. To find the maximum value function, The function is for finding the minimum value; To round up; The target difference between two adjacent vertical pressure levels is set by estimating the internal friction angle of the infill material based on regional geological experience, and then setting different target differences based on different estimated internal friction angles of the infill material. The method for setting k is based on standard experimental statistical requirements to ensure data reliability and error assessment, and a preset shear rate is used. The setting method is as follows: it is set according to the permeability coefficient and height of the filling material, such as... ,in, The permeability coefficient is obtained primarily through variable head permeability tests. The specific measurement method is a conventional technique and will not be described in detail here. The height of the filling material; the method for determining shear loss of the sample filling material is as follows: if any one of the criteria is met, it is determined to be shear loss, specifically:

[0117] Criterion 1: Continuously measure the force gauge readings. When the absolute value of the change between two consecutive readings in N1 consecutive readings is less than or equal to the first change threshold, and the difference between the maximum and minimum values ​​in the corresponding N1 consecutive readings is less than or equal to the second change threshold, it is determined to be shear loss. The method for setting N1 is: based on the frequency of use and the random noise when the filling material is sheared. The method for setting the first change threshold and the second change threshold is: the first change threshold can be set according to the minimum scale division value of the force gauge dial indicator, and the second change threshold is set by adding the actual measurement error of the force gauge to the first change threshold.

[0118] Criterion 2: Record the maximum reading during the shearing process. When the reading drops from the maximum reading by a factor greater than or equal to the magnitude threshold, it is considered shear loss. The magnitude threshold is set based on the maximum reading and engineering experience. Generally, a corresponding disturbance is set based on the maximum reading. Combined with engineering experience, when the shear strength decreases by more than 5%, that is, the maximum reading multiplied by 5%, it indicates that the material has failed.

[0119] Criterion 3: Install a displacement sensor to measure the shear displacement. When the shear displacement reaches the displacement threshold, it is determined to be shear loss even if neither Criterion 1 nor Criterion 2 is triggered. The displacement threshold is set as the maximum allowable shear displacement of the direct shear instrument when the horizontal shear force is applied.

[0120] Filler cohesion when activation sensitivity type is A The setting method is as follows: in the setting method of the cohesion of the filling material corresponding to the internal friction angle of the filling material, the intercept of the straight line of the maximum shear stress with respect to the vertical pressure obtained by fitting is on the vertical axis; the rock mass integrity coefficient when the sensitivity type is activated is A. The calculation method includes: if core quality data exists, the rock mass integrity coefficient is obtained by calculating the ratio of the core quality index value to the preset value; if core quality data is lacking, the rock mass integrity coefficient is obtained by calculating the square of the ratio of the P-wave velocity to the P-wave velocity of the intact rock mass; wherein, the P-wave velocity of the intact rock mass is determined according to relevant standards and in combination with lithology type; the preset value is the core quality index value of the intact rock mass, which is set according to the core quality index benchmark value corresponding to the intact rock mass in the rock mass integrity evaluation standard.

[0121] For type B, which corresponds to the moderately sensitive type, the instability mode is shear failure of cohesive soil along the structural plane. Based on stress wave theory and elastic wave theory, the critical vibration threshold for activating the sensitivity type B is... ,in, To activate the cohesive force of the filler when the sensitivity type is B; The thickness of the waterproof roof when the sensitivity type is B; To activate the internal friction angle of the filler when the sensitivity type is B; To activate the rock mass integrity coefficient when the sensitivity type is B; The natural density of the filler when the activation sensitivity type is B; The parameter acquisition method for calculating the critical vibration threshold when the activation sensitivity type is B is the same as the acquisition method for the cohesion of the filling material when the activation sensitivity type is A.

[0122] For Class C, which corresponds to the low-sensitivity type, the instability mode is crack propagation caused by vibration, but there is no risk of instability of the filling medium. The main control point is to avoid the opening of the rock mass structure surface. The critical vibration threshold for Class C can be set according to relevant industry standards, such as taking the allowable vibration velocity limit corresponding to the stability of the surrounding rock of the tunnel, the blasting type and the protected object in the "Blasting Safety Regulations". When there are multiple allowable vibration velocity limits, the smaller value is taken as the critical vibration threshold for Class C.

[0123] The blasting design parameters of each blasting point in the blasting operation to be carried out are collected. Each blasting point in the blasting operation to be carried out is screened to obtain effective blasting points. Based on the blasting design parameters of the effective blasting points, the predicted peak velocity generated by the blasting operation at the target measuring point is predicted. Among them, the blasting design parameters include the single-stage detonation charge, the location of the blasting point, and the detonation delay time corresponding to each blasting point. The blasting operation to be carried out is the blasting operation that will be carried out at the working face or in the tunnel construction area. The blasting point is the location of the blasting source corresponding to each blast hole or each detonation segment in the blasting operation to be carried out; the target measuring point is the location of the concealed karst water body, filling medium, or karst anomaly; the spatial straight-line distance between the blasting point and the target measuring point, i.e., the blast center distance, is used to characterize the distance attenuation of the blasting vibration as it propagates from the blasting source to the target measuring point; a blasting point whose distance from the target measuring point is less than or equal to the preset influence distance refers to a blasting point whose blasting vibration may propagate to the target measuring point and cause disturbance to the concealed karst water body or filling medium at the target measuring point; the preset influence distance is obtained by back-calculation according to the Sadovsky formula, that is, using the maximum single-stage detonation charge, geological condition correction coefficient, and attenuation index of the blasting operation to be carried out as inputs, and using the minimum value among the critical vibration thresholds corresponding to each activation sensitivity type as the control peak velocity, the blast center distance is obtained by back-calculation, and the corresponding blast center distance is used as the preset influence distance.

[0124] Based on the disturbance source parameters, the blasting points in the blasting operation to be carried out are screened to obtain the effective blasting points that meet the disturbance influence conditions of the target measuring point. Based on the disturbance source parameters corresponding to the effective blasting points, the predicted disturbance intensity of the blasting operation to be carried out at the target measuring point is predicted.

[0125] In a specific implementation, refer to Figure 2 Methods for obtaining the predicted disturbance intensity at the target measurement point include:

[0126] Obtain the blasting design parameters for each blasting point in the blasting operation to be carried out;

[0127] For any target measuring point, select the blasting points in the blasting operation to be carried out whose distance from the target measuring point is less than or equal to the preset influence distance, and use them as the effective blasting points corresponding to the target measuring point;

[0128] Calculate the spatial straight-line distance between each effective blasting point and the target measuring point to obtain the blast center distance of each effective blasting point relative to the target measuring point;

[0129] Obtain the geological condition correction coefficient corresponding to the target measuring point, and match the geological condition correction coefficient from the preset coefficient mapping table based on the rock mass integrity coefficient of the rock mass between the working face and the target measuring point.

[0130] The methods for calculating the rock mass integrity factor include:

[0131] If core quality data is available, the rock mass integrity coefficient is obtained by calculating the ratio of the core quality index value to the preset value; if core quality data is lacking, the rock mass integrity coefficient is obtained by calculating the square of the ratio of the P-wave velocity to the P-wave velocity of the intact rock mass; wherein, the P-wave velocity of the intact rock mass is determined according to relevant standards and in combination with the lithology type.

[0132] The method for establishing the preset coefficient mapping table is as follows: Based on the rock mass integrity coefficient, the rock mass is divided into four levels, each level corresponding to a geological condition correction coefficient interval, and the midpoint of the interval is taken as the mapping value; specifically:

[0133] When the rock mass integrity coefficient is not less than the first coefficient threshold, it is marked as an intact rock mass, corresponding to the first geological condition correction coefficient;

[0134] When the rock mass integrity coefficient is less than the first coefficient threshold but not less than the second coefficient threshold, it is marked as a relatively intact rock mass, corresponding to the second geological condition correction coefficient.

[0135] When the rock mass integrity coefficient is less than the second coefficient threshold but not less than the third coefficient threshold, it is marked as a fractured rock mass, corresponding to the third geological condition correction coefficient.

[0136] When the rock mass integrity coefficient is less than the third coefficient threshold, it is marked as an extremely fractured rock mass, corresponding to the fourth geological condition correction coefficient. The first, second, and third coefficient thresholds are used to map the rock mass integrity coefficient to different geological condition correction coefficients. If core quality data exists, the rock mass integrity coefficient is calculated as the ratio of the core quality index value to the benchmark value of the intact core quality index. If core quality data is missing, the rock mass integrity coefficient is calculated as the square of the ratio of the P-wave velocity at the target measuring point to the P-wave velocity of the intact rock mass. The first, second, and third coefficient thresholds are determined based on the geological condition correction calibration sample set. Each sample in the geological condition correction calibration sample set includes the rock mass integrity coefficient, single-stage detonation charge, detonation center distance, and measured peak velocity.

[0137] For each sample in the geological condition correction calibration sample set, the equivalent geological condition correction coefficient corresponding to that sample is calculated using the blasting vibration propagation formula. In one embodiment, the equivalent geological condition correction coefficient... The calculation formula is: ;in, This represents the measured peak velocity of the corresponding sample. This represents the amount of explosive charge for a single stage of detonation corresponding to this sample. This represents the distance from the center of the explosion corresponding to this sample. The decay index is obtained by regression based on the same calibrated sample set.

[0138] The geological condition correction calibration sample set is sorted from largest to smallest according to the rock mass integrity coefficient. Within the sorted sample sequence, three candidate coefficient boundary points are selected to divide the sample sequence into four consecutive sample groups. For each candidate division method, the intra-group dispersion of the equivalent geological condition correction coefficients within the four consecutive sample groups is calculated, and the sum of the intra-group dispersions of the four consecutive sample groups is taken as the division error of that candidate division method.

[0139] Within-group dispersion is calculated as the sum of the squared differences between the equivalent geological condition correction coefficients of each sample in the group and the average of the equivalent geological condition correction coefficients of the group.

[0140] The candidate coefficient boundary point with the smallest division error is determined as the target division method. The first coefficient threshold is the arithmetic mean of the rock mass integrity coefficients of adjacent samples on both sides of the boundary between the first and second groups; the second coefficient threshold is the arithmetic mean of the rock mass integrity coefficients of adjacent samples on both sides of the boundary between the second and third groups; and the third coefficient threshold is the arithmetic mean of the rock mass integrity coefficients of adjacent samples on both sides of the boundary between the third and fourth groups.

[0141] When there are two or more candidate partitioning methods with the same partitioning error, the candidate partitioning method that meets the following conditions shall be selected: the smaller the rock mass integrity coefficient, the larger the average value of the equivalent geological condition correction coefficient of the corresponding sample group; if there are still two or more candidate partitioning methods, the candidate partitioning method with the smallest difference in the number of samples among the four consecutive sample groups shall be selected.

[0142] The method for setting the first, second, third, and fourth geological condition correction coefficients is as follows: the range is based on the statistical analysis of a large amount of blasting vibration monitoring data. Regression is performed on the blasting vibration monitoring data of similar karst tunnels to obtain the distribution of geological condition correction coefficients corresponding to each rock mass integrity interval. The median or arithmetic mean of the distribution of geological condition correction coefficients in the corresponding interval is used as the mapping value.

[0143] Based on the Sadovsky formula, and combined with geological condition correction factors, single-stage detonation charge, and detonation center distance, the segmented predicted peak velocity generated at the target measuring point by each effective blasting point is calculated. In one feasible implementation, the segmented predicted peak velocity can be calculated using the following formula: In, among them, For the first The effective blasting point is at the first The segmented predicted peak velocity generated at each target measurement point; For the first Geological condition correction coefficients for each target measuring point; For the first The attenuation index corresponding to each effective blasting point was obtained through regression analysis of a large amount of blasting vibration monitoring data. For the first The amount of explosive charge per effective blasting point; For the first From the first effective blasting point to the... The distance between the explosion center and the target measurement point.

[0144] When there are multiple target measurement points, each target measurement point is treated as an independent risk assessment object, and the corresponding predicted peak velocity is calculated one by one. The predicted peak velocities between different target measurement points are not superimposed.

[0145] For any target measuring point, if the time difference between the detonation delay of each effective blasting point is greater than the preset vibration duration threshold, then the maximum value of the predicted peak velocity of each segment is taken as the predicted peak velocity at the target measuring point, and the predicted peak velocity is taken as the predicted disturbance intensity. The preset vibration duration threshold is set by determining the duration of a single segment of detonation vibration signal from the start of blasting to the decay to the background vibration level in the on-site blasting vibration monitoring waveform.

[0146] If the time difference between the initiation delay of at least two effective blasting points is less than or equal to the preset vibration duration threshold, the segmented predicted peak velocities corresponding to the at least two effective blasting points are conservatively synthesized to obtain the composite predicted peak velocity. The maximum value among the composite predicted peak velocity and the remaining segmented predicted peak velocities is taken as the predicted peak velocity at the target measuring point. The conservative synthesis method is as follows: the segmented predicted peak velocities corresponding to each effective blasting point with an initiation delay time difference less than or equal to the preset vibration duration threshold are added together to obtain the composite predicted peak velocity.

[0147] By comparing and analyzing the predicted disturbance intensity and the critical disturbance threshold, the activation risk coefficient of the hidden karst water body at the target measuring point is obtained.

[0148] In specific implementation methods, the activation risk coefficient of concealed karst water bodies at the target measuring point includes:

[0149] The activation risk coefficient is calculated based on the critical vibration threshold and the predicted peak velocity corresponding to the target measuring point. That is, the predicted peak velocity is divided by the critical vibration threshold to obtain the activation risk coefficient.

[0150] When the activation risk coefficient is greater than the preset risk threshold, it indicates that the hidden karst water body at the target measuring point is at risk of being activated by blasting vibration, and microseismic precursor identification needs to be carried out during the blasting operation; when the activation risk coefficient is less than or equal to the preset risk threshold, it indicates that the hidden karst water body at the target measuring point is under routine monitoring; the preset risk threshold is a safety margin coefficient less than or equal to 1, which is determined based on the deviation distribution between the predicted peak velocity and the measured peak velocity in similar tunnel blasting monitoring data; when the prediction deviation is large, the preset risk threshold is reduced.

[0151] If the activation risk coefficient is greater than the preset risk threshold, microseismic activity characteristics are generated based on the blasting vibration response data, and precursor characteristic indicators are generated based on the activation sensitivity type and microseismic activity characteristics.

[0152] In a specific implementation, the method for generating precursor characteristic indicators includes:

[0153] When the activation risk coefficient is greater than the preset risk threshold, blasting vibration response data corresponding to the target measuring point is collected during the blasting operation. Microseismic activity characteristics are obtained based on the analysis of the blasting vibration response data. Specifically, microseismic data is continuously collected from the first preset time before the blasting operation until the second preset time after the blasting. The microseismic data includes the occurrence time, energy density, and spectral distribution of microseismic events. The proportion of high-frequency microseismic signals is obtained based on the spectral distribution. The method for setting the first and second preset times is determined based on the coverage requirements of the instability response process of the filling medium before and after the blasting.

[0154] When the microseismic frequency anomaly rate is not lower than the frequency anomaly threshold of the corresponding active sensitive type at the target measuring point, the precursor features of the microseismic frequency anomaly are triggered; otherwise, the precursor features of the microseismic frequency anomaly are not triggered. The method for setting the frequency anomaly threshold is as follows: blasting tests with different amounts of explosives can be carried out in sections with known active sensitive types, the response of the filling medium is monitored simultaneously, the frequency anomaly rate corresponding to the start of micro-fractures in the filling medium is recorded, and the corresponding frequency anomaly rate is taken as the frequency anomaly threshold.

[0155] When the energy density anomaly rate is not lower than the density anomaly threshold of the corresponding activation sensitive type of the target measuring point, the precursor features of energy density anomaly are triggered; otherwise, the precursor features of energy density anomaly are not triggered. The method for setting the density anomaly threshold is as follows: In the blasting test in the early stage of construction, the changes in water level and seepage flow in the borehole are monitored simultaneously. When the seepage flow first shows a continuous increase and the increase exceeds 5 times the initial seepage value within 30 minutes, the corresponding microseismic energy density value is recorded. 80% of the corresponding microseismic energy density value is used as the lower limit of density anomaly warning to avoid the uncertainty of adopting a single test value.

[0156] When the proportion of high-frequency signals is not lower than the signal proportion threshold of the corresponding activation sensitive type at the target measuring point, the precursor features of the high-frequency signal are triggered; otherwise, the precursor features of the high-frequency signal are not triggered. The method for setting the signal proportion threshold is as follows: a blasting test can be carried out in a section with a known filling medium type, and the response of the filling medium, such as changes in seepage flow, can be monitored simultaneously. The signal proportion value corresponding to when the filling medium begins to show abnormal seepage is recorded, and the corresponding signal proportion value is taken as the signal proportion threshold.

[0157] The risk identification results of sudden water inrush corresponding to the target measuring point are generated based on the activation risk coefficient and precursor characteristic indicators.

[0158] In a specific implementation, refer to Figure 3 Methods for generating the risk identification results of sudden water inrush corresponding to the target measuring point include:

[0159] Based on the trigger status of the precursor features, a real-time early warning signal is generated according to the following rules: Specifically, if the number of triggered precursor features is not less than 2, a Level 1 early warning signal is generated, prompting the blasting operation to be stopped, personnel to be evacuated to a safe area, and the emergency plan to be activated.

[0160] If the number of precursor features that are triggered is 1, a level 2 warning signal is generated, prompting a reduction in the amount of subsequent explosives and an increase in the monitoring frequency;

[0161] If the number of precursor features triggered is 0, a level 3 warning signal is generated, prompting the user to maintain encrypted monitoring and emergency equipment standby status.

[0162] If the activation risk coefficient corresponding to the target measuring point is less than or equal to the preset risk threshold, microseismic precursor identification will not be triggered, and a regular monitoring signal will be generated.

[0163] Based on real-time early warning signals, when a first-level or second-level early warning signal is generated, the blasting parameters are optimized and adjusted to control the blasting vibration below the critical vibration threshold.

[0164] In specific implementation methods, the methods for optimizing and adjusting blasting parameters include:

[0165] Based on the Sadovsky formula, combined with geological condition correction coefficients, the critical vibration threshold corresponding to the target measuring point, and the distance from the detonation center, the maximum allowable single-stage charge is calculated. In one feasible implementation, the maximum allowable single-stage charge can be calculated using the following formula: ;in, For the first The first explosion point is relative to the first The maximum allowable single-stage drug dosage for each target measurement point; For the first The critical vibration threshold corresponding to each target measuring point; when the same blasting point corresponds to multiple target measuring points, the minimum value among the maximum allowable single-stage charge calculated for each target measuring point is taken as the optimized maximum single-stage charge for the blasting point, so that the optimized blasting parameters simultaneously meet the critical vibration control requirements of multiple target measuring points; the single-stage charge for the corresponding blasting point in the blasting operation to be carried out is controlled within the maximum single-stage charge to obtain the adjusted maximum single-stage charge.

[0166] Example 2

[0167] Data acquisition module: acquires multi-source input data corresponding to the target measuring point under the blasting disturbance scenario of karst tunnel. The multi-source input data includes the structural state data of the target measuring point, the disturbance source parameters of the blasting operation to be carried out, and the blasting vibration response data.

[0168] Sensitive Classification Module: Based on structural state data analysis, the activation sensitivity type corresponding to the filling medium at the target measurement point is obtained, and the corresponding critical disturbance threshold is matched from the preset threshold mapping table according to the activation sensitivity type;

[0169] Disturbance prediction module: Based on the disturbance source parameters, the blasting points in the blasting operation to be carried out are screened to obtain the effective blasting points that meet the disturbance influence conditions of the target measuring point, and the predicted disturbance intensity of the blasting operation to be carried out at the target measuring point is predicted based on the disturbance source parameters corresponding to the effective blasting points.

[0170] Comparison and Judgment Module: Based on the predicted disturbance intensity and the critical disturbance threshold, a comparative analysis is performed to obtain the activation risk coefficient of the hidden karst water body at the target measuring point;

[0171] Precursor analysis module: If the activation risk coefficient is greater than the preset risk threshold, microseismic activity characteristics are generated based on the blasting vibration response data, and precursor characteristic indicators are generated based on the activation sensitivity type and microseismic activity characteristics.

[0172] Analysis and early warning module: Generates the risk identification results of sudden water inrush corresponding to the target measuring point based on the activation risk coefficient and precursor characteristic indicators.

[0173] Finally: The above description is only a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A method for multi-source data risk identification in blasting disturbance scenarios, characterized in that, include: Acquire multi-source input data corresponding to the target measuring point under the blasting disturbance scenario of karst tunnel. The multi-source input data includes the structural state data of the target measuring point, the disturbance source parameters of the blasting operation to be carried out, and the blasting vibration response data. Based on the type of filling medium and the thickness of the impermeable roof in the structural state data, and combined with the preset type determination conditions, the activation sensitivity type corresponding to the filling medium at the target measuring point is determined by analysis: Based on the type of filling medium and the thickness of the impermeable roof in the structural state data of the target measuring point, if the filling medium type is sandy soil, or the thickness of the impermeable roof is less than the first thickness threshold, then the activation sensitivity type corresponding to the filling medium at the target measuring point is determined to be type A. If the filling medium is cohesive soil, or the thickness of the waterproof top plate is greater than or equal to the first thickness threshold and less than or equal to the second thickness threshold, then the activation sensitivity type corresponding to the filling medium at the target measuring point is determined to be Class B. If the filling medium type is no filling, or the thickness of the impermeable roof is greater than the second thickness threshold, then the activation sensitivity type corresponding to the filling medium at the target measuring point is determined to be Class C; among them, the activation sensitivity of Class A, Class B and Class C decreases in that order; when the determination result corresponding to the filling medium type and the thickness of the impermeable roof is inconsistent, the activation sensitivity type with the higher risk level is taken as the activation sensitivity type corresponding to the filling medium at the target measuring point; according to the activation sensitivity type, the corresponding critical disturbance threshold is matched from the preset threshold mapping table; Based on the disturbance source parameters, the blasting points in the blasting operation to be carried out are screened to obtain the effective blasting points that meet the disturbance influence conditions of the target measuring point. Based on the disturbance source parameters corresponding to the effective blasting points, the predicted disturbance intensity of the blasting operation to be carried out at the target measuring point is predicted. Based on the comparative analysis of the predicted disturbance intensity and the critical disturbance threshold, the activation risk coefficient of the hidden karst water body at the target measuring point is obtained: the ratio of the predicted peak velocity corresponding to the target measuring point to the critical vibration threshold corresponding to the target measuring point is calculated to obtain the activation risk coefficient of the hidden karst water body at the target measuring point. If the activation risk coefficient is greater than the preset risk threshold, then microseismic activity characteristics are generated based on the blasting vibration response data, and precursor characteristic indicators are generated based on the activation sensitivity type and microseismic activity characteristics: Based on the activation sensitivity type corresponding to the target measuring point, match the corresponding frequency anomaly threshold, density anomaly threshold, and signal proportion threshold; based on the microseismic data collected within the first preset time before blasting, obtain the background value of microseismic event frequency and the background value of microseismic energy density; based on the microseismic data collected within the second preset time after blasting, obtain the microseismic event frequency, microseismic energy density, and proportion of high-frequency microseismic signals; calculate the microseismic frequency anomaly rate based on the background value of microseismic event frequency and the microseismic event frequency, and calculate the energy density anomaly rate based on the background value of microseismic energy density and the microseismic energy density; use the microseismic frequency anomaly rate, energy density anomaly rate, and proportion of high-frequency microseismic signals as precursor characteristic indicators; The risk identification results of sudden water inrush corresponding to the target measuring point are generated based on the activation risk coefficient and precursor characteristic indicators.

2. The method for multi-source data risk identification in blasting disturbance scenarios according to claim 1, characterized in that, Methods for obtaining the activation sensitivity type corresponding to the filling medium at the target measurement point also include: If the type of filling medium and the thickness of the waterproof roof do not meet the preset type determination conditions, the core quality index value is obtained by analyzing the drilling data, and ground-penetrating radar data and tunnel seismic wave advance prediction data are collected. Ground-penetrating radar characteristics and tunnel seismic wave advance prediction characteristics are obtained by analyzing the ground-penetrating radar data and tunnel seismic wave advance prediction data. The activation sensitive type corresponding to the filling medium at the target measuring point is obtained by jointly determining the core quality index value, ground-penetrating radar characteristics and tunnel seismic wave advance prediction characteristics.

3. The method for multi-source data risk identification in blasting disturbance scenarios according to claim 1, characterized in that, Methods for obtaining the predicted disturbance intensity at the target measuring point during a blasting operation include: Collect blasting design parameters for the blasting points to be blasted, including the location of the blasting point, the amount of explosive charge per stage, and the detonation delay time; screen blasting points whose distance from the blasting point to be blasted to the target measuring point is less than or equal to the preset influence distance, and use them as valid blasting points corresponding to the target measuring point; Calculate the spatial straight-line distance between each effective blasting point and the target measuring point to obtain the blast center distance of each effective blasting point relative to the target measuring point; Based on the single-stage detonation charge, detonation center distance and geological condition correction coefficient corresponding to each effective blasting point, the segmented predicted peak velocity generated at the target measuring point at each effective blasting point is calculated respectively. Based on the analysis of the predicted peak velocity of each segment, the predicted peak velocity generated at the target measuring point by the blasting operation to be carried out is obtained, and the predicted peak velocity is used as the predicted disturbance intensity.

4. The method for multi-source data risk identification in blasting disturbance scenarios according to claim 3, characterized in that, Methods for obtaining the predicted peak velocity at the target measuring point based on the analysis of predicted peak velocities for each segment include: If it is determined that the initiation delay time difference between each effective blasting point is greater than the preset vibration duration threshold, then the maximum value among the predicted peak velocities of each segment is taken as the predicted peak velocity at the target measuring point. If the time difference between the detonation delay between at least two effective detonation points is less than or equal to a preset vibration duration threshold, then the segmented predicted peak velocities corresponding to the at least two effective detonation points are conservatively synthesized to obtain a synthesized predicted peak velocity, and the maximum value among the synthesized predicted peak velocity and the remaining segmented predicted peak velocities is taken as the predicted peak velocity at the target measuring point.

5. The method for multi-source data risk identification in blasting disturbance scenarios according to claim 1, characterized in that, Methods for obtaining the risk identification results of sudden water inrush corresponding to the target measuring point include: When the microseismic frequency anomaly rate is not lower than the frequency anomaly threshold of the corresponding activation sensitive type of the target measuring point, the precursor features of the microseismic frequency anomaly are triggered. When the energy density anomaly rate is not lower than the density anomaly threshold of the activation sensitive type corresponding to the target measurement point, the precursor features of energy density anomaly are triggered. When the proportion of high-frequency microseismic signals is not lower than the threshold of the proportion of signals of the corresponding active sensitive type at the target measuring point, the precursor features of the high-frequency signal are triggered. The risk identification result of sudden water inrush corresponding to the target measuring point is generated based on the number of precursor features that are triggered.

6. The method for multi-source data risk identification in blasting disturbance scenarios according to claim 5, characterized in that, Methods for generating the risk identification results of sudden water inrush at target measuring points based on the number of triggered precursor features include: If the number of precursor features that are triggered is not less than 2, a Level 1 warning signal is generated; If the number of precursor features that are triggered is 1, a level 2 warning signal is generated; If the number of precursor features that are triggered is 0, a Level 3 warning signal is generated.

7. A multi-source data risk identification system for blasting disturbance scenarios, implementing the multi-source data risk identification method for blasting disturbance scenarios as described in any one of claims 1-6, characterized in that, include: Data acquisition module: acquires multi-source input data corresponding to the target measuring point under the blasting disturbance scenario of karst tunnel. The multi-source input data includes the structural state data of the target measuring point, the disturbance source parameters of the blasting operation to be carried out, and the blasting vibration response data. Sensitive Classification Module: Based on structural state data analysis, the activation sensitivity type corresponding to the filling medium at the target measurement point is obtained, and the corresponding critical disturbance threshold is matched from the preset threshold mapping table according to the activation sensitivity type; Disturbance prediction module: Based on the disturbance source parameters, the blasting points in the blasting operation to be carried out are screened to obtain the effective blasting points that meet the disturbance influence conditions of the target measuring point, and the predicted disturbance intensity of the blasting operation to be carried out at the target measuring point is predicted based on the disturbance source parameters corresponding to the effective blasting points. Comparison and Judgment Module: Based on the predicted disturbance intensity and the critical disturbance threshold, a comparative analysis is performed to obtain the activation risk coefficient of the hidden karst water body at the target measuring point; Precursor analysis module: If the activation risk coefficient is greater than the preset risk threshold, microseismic activity characteristics are generated based on the blasting vibration response data, and precursor characteristic indicators are generated based on the activation sensitivity type and microseismic activity characteristics. Analysis and early warning module: Generates the risk identification results of sudden water inrush corresponding to the target measuring point based on the activation risk coefficient and precursor characteristic indicators.

Images