Tunnel surrounding rock ground stress intelligent analysis method and system based on while-drilling parameters

Abstract

The application relates to the technical field of tunnel engineering geology construction, and discloses a tunnel surrounding rock ground stress intelligent analysis method and system based on a while-drilling parameter, which comprises the following steps: through multi-source data fusion and mixed model cooperation, dynamic coordinate system construction, mechanical specific energy calculation, Stacking model classification, wave velocity ellipsoid modeling, dynamic weight fusion and difference-driven re-calibration are realized, and three-dimensional stress field high-precision reconstruction and construction parameter optimization are completed; the system comprises a multi-source sensing module, a dynamic modeling module, a mixed analysis module, a self-adaptive fusion module, a three-dimensional reconstruction module and a closed-loop control module. Through multi-source data fusion and mixed model cooperation, combined with a dynamic weight distribution and a local re-calibration mechanism, high-precision ground stress analysis is realized. The wave velocity ellipsoid physical model and the Stacking integrated learning model are complementarily used, combined with octree space reconstruction and closed-loop control, the adaptability to complex geology is improved, and the construction efficiency and system reliability are optimized.

Description

Technical Field

[0001] This invention relates to the field of tunnel engineering geological construction technology, specifically to a method and system for intelligent analysis of tunnel surrounding rock stress based on drilling parameters. Background Technology

[0002] In tunnel engineering, accurate analysis of the surrounding rock stress field is crucial for ensuring construction safety and optimizing support design. Traditional in-situ stress testing relies heavily on point-based measurement methods such as borehole stress relief and hydraulic fracturing, which have limitations including long testing cycles, high costs, and the inability to acquire continuous stress field data in real time. While measurement-while-drilling (MWD) technology enables dynamic acquisition of drilling process parameters, existing methods often employ single physical or statistical models for stress analysis, making it difficult to effectively integrate multi-source heterogeneous data characteristics and prone to model mismatch issues under complex geological conditions.

[0003] Current geostress analysis systems based on drilling parameters generally suffer from the following technical bottlenecks: Physically driven models rely excessively on theoretical assumptions, resulting in insufficient adaptability to complex working conditions such as rock mass anisotropy and fracture development; data-driven models lack interpretability, easily leading to erroneous inferences in areas not covered by training data; multi-model fusion strategies often employ fixed weight allocation, failing to dynamically respond to abrupt changes in geological conditions; and spatial continuity constraints are ignored during 3D stress field reconstruction, leading to missed detection of local stress concentration phenomena. Furthermore, existing systems lack closed-loop control interfaces with construction machinery, making it difficult to provide real-time decision support for intelligent construction.

[0004] Therefore, this invention proposes an intelligent analysis method and system for tunnel surrounding rock stress based on drilling parameters to address the shortcomings of existing technologies. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides an intelligent analysis method and system for tunnel surrounding rock stress based on drilling parameters, which solves the problems of insufficient multi-source data fusion, poor model adaptability, and lack of real-time closed-loop control in existing tunnel surrounding rock stress analysis methods.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a smart analysis method for tunnel surrounding rock stress based on drilling parameters, the method comprising the following steps:

[0007] S1. Real-time acquisition of drilling parameters, triaxial wave velocity data and rock mass quality parameters through sensor array;

[0008] S2. Establish a dynamic local coordinate system based on the propulsion speed vector direction in the drilling parameters, and calculate the mechanical specific energy according to the impact pressure, rotation torque and rock breaking volume in the drilling parameters.

[0009] S3. Input the mechanical specific energy and the rock mass quality parameters into the Stacking ensemble learning model to classify the geostress state and output the rock strength stress ratio classification result.

[0010] S4. Based on the geostress state classification results and the triaxial wave velocity data, analyze the three-dimensional wave velocity distribution and construct a wave velocity ellipsoid physical model by fitting the least squares method.

[0011] S5. By integrating the analytical results of the wave velocity ellipsoid physical model with the analytical results of the Stacking ensemble learning model, the principal stress direction is calculated using a dynamic weight allocation strategy based on rock mass integrity.

[0012] S6. When the relative difference between the analytical results of the wave velocity ellipsoid physical model and the analytical results of the Stacking ensemble learning model exceeds a set threshold, a local recalibration mechanism is triggered and the three-dimensional stress field model is updated.

[0013] Preferably, in step S1:

[0014] The sensor group includes a drilling parameter sensor, a triaxial ultrasonic probe, and a ground radar installed on the intelligent rock drilling rig.

[0015] The drilling parameters include the impact pressure P, which characterizes the working state of the rock drilling machinery. imp , slewing torque τ, propulsion speed vector v drill Drilling diameter D h and piston rod structural parameters D r ;

[0016] The triaxial wave velocity data is acquired through an ultrasonic array probe arranged around the borehole, including the x-axis wave velocity V in an orthogonal coordinate system. x y-axis wave velocity V y z-axis wave velocity V z , where the z-axis is aligned with the drilling direction;

[0017] The rock mass quality parameter is the rock quality index RQD obtained through ground-penetrating radar scanning, and its calculation satisfies:

[0018]

[0019] Among them, l intact L represents the length of a complete rock mass segment within the scanned area. scan This represents the total scan length.

[0020] The drilling parameters, triaxial wave velocity data, and rock mass quality parameters are stored in a buffer queue in a timestamp-aligned manner.

[0021] Preferably, the method for establishing the dynamic local coordinate system in step S2 is as follows:

[0022] The starting point of the drilling operation is taken as the spatial reference origin. The propulsion velocity vector v drill The unit direction vector is defined as the positive z-axis direction. The geometric surface S of the face of the machine is obtained by laser scanning, and its normal vector is calculated. As the y-axis direction, the x-axis direction is generated by the y×z orthogonal direction determined by the right-hand rule;

[0023] The calculation of the mechanical specific energy is based on the combined effect of impact rock-breaking energy and rotary cutting energy, and its expression is:

[0024]

[0025] in, L represents the rock-breaking volume in a single drilling cycle. pen Drilling depth; L is the hydraulic impact piston stroke; P imp Peak impact pressure; D h D is the borehole diameter; r τ is the diameter of the impact piston rod; τ is the output torque of the rotary motor; v pen This represents the real-time drilling speed.

[0026] Preferably, in the mechanical specific energy calculation:

[0027] Impact rock-breaking energy term Reflects the working efficiency of the hydraulic impact system;

[0028] Rotary cutting energy term Characterizing the rotational shear effect of the cutting tool in rock breaking;

[0029] Rock breaking volume V rock The calculation introduces the borehole diameter D h The square term is used to reflect the effect of the pore wall area on the energy density, and the unit of the mechanical specific energy SE is J / m². 3 .

[0030] Preferably, in step S3:

[0031] The Stacking ensemble learning model adopts a two-level architecture. The first-level base model group includes support vector machine, K-nearest neighbor algorithm and random forest. The second-level meta-model uses linear regression algorithm to weight and fuse the outputs of the base models.

[0032] The rock strength-stress ratio classification result is calculated using the following formula:

[0033]

[0034] Where n is the total number of geostress state categories; p i λ represents the probability of the i-th type of geostress state determined by the model.i These are the feature coefficients for the corresponding categories;

[0035] The characteristic coefficients are pre-calibrated based on the mapping relationship between different rock mass strengths and stress states in historical borehole data, and the input feature vector includes normalized mechanical specific energy. Rock quality indicators and the interaction term of the three-axis wave velocity parameters

[0036] Preferably, in step S4:

[0037] The three-dimensional wave velocity analysis is based on the coupling relationship between mechanical specific energy and the anisotropic characteristics of rock mass, and a wave velocity prediction model is established:

[0038] V i =α i SE 2 +β i SE+γ i ,(i∈{x,y,z});

[0039] Where, α i β i γ i The correlation coefficient reflects the elastic modulus, Poisson's ratio, and degree of fracture development of the rock mass;

[0040] The physical model of the wave velocity ellipsoid is constructed by optimizing the ellipsoid parameters using the least squares method, with the objective function being:

[0041]

[0042] Where A, B, and C represent the semi-axis lengths of the ellipsoid in the x, y, and z axes, respectively, and their physical meaning corresponds to the wave velocity distribution characteristics in the direction of maximum principal stress; the regularization term λ(A 2 +B 2 +C 2 Control the complexity of the model to prevent overfitting.

[0043] Preferably, in step S5:

[0044] The dynamic weight allocation strategy adjusts the fusion weights of the physical model and the AI ​​model using two factors: rock mass integrity and stress gradient. The calculation formula is as follows:

[0045]

[0046] Among them, w p η is the weighting coefficient of the wave velocity ellipsoid physical model; RQD0 is the preset rock mass integrity benchmark threshold; k is the integrity adjustment factor, which controls the sensitivity of the weights to changes in RQD; η is the stress gradient adjustment factor, which controls the strength of the weights' response to stress abrupt changes. The stress field gradient magnitude is obtained by calculating the partial derivatives of the three-dimensional stress tensor.

[0047] The final principal stress direction analytical value is calculated by the following formula:

[0048] σ final =w p σ phys +(1-w p )σ AI ;

[0049] Where, σ phys This represents the analytical results of the wave velocity ellipsoid physical model; σ AI This represents the predicted output of the Stacking ensemble learning model.

[0050] Preferably, in step S6:

[0051] The determination of the relative difference adopts the normalized error assessment method, and its calculation formula is:

[0052]

[0053] When δ > ε, a local recalibration mechanism is triggered, where ε is a preset difference threshold.

[0054] The recalibration process takes place in the abnormal subspace partitioned by the octree. Internally, a parameter optimization algorithm based on historical data similarity is used to update the wave velocity ellipsoid model coefficients α. i ,β i γ i And the weight parameters of the Stacking ensemble learning model.

[0055] Preferably, the octree subspace The division is dynamically generated based on the spatial topological relationship between lidar point cloud data and real-time drilling trajectory, and its boundary conditions satisfy:

[0056]

[0057] in, Boundary values ​​are determined using a point cloud density clustering algorithm;

[0058] The three-dimensional stress field model is updated using the B-spline interpolation algorithm, the expression of which is:

[0059] σ(x,y,z)=∑ i,j,k B i,d (x)B j,d (y)B k,d (z)σ ijk ;

[0060] Where σ(x,y,z) is the three-dimensional stress field interpolation function; B i,d σ is a B-spline basis function of order d; ijk The stress value at the interpolation control point is given, and the interpolation order d≥3 is used to ensure the continuity of the stress field.

[0061] This invention also provides an intelligent analysis system for tunnel surrounding rock stress based on drilling parameters, the system comprising the following modules:

[0062] The multi-source sensing module integrates a drilling-while-drilling sensor array, a three-axis ultrasonic probe, and a ground-penetrating radar, enabling the synchronous acquisition of drilling parameters, wave velocity data, and rock mass quality parameters.

[0063] The dynamic modeling module includes a local coordinate system generation unit and a mechanical energy calculation unit, which constructs drilling space benchmarks and energy characteristics in real time.

[0064] The hybrid analysis module deploys a Stacking ensemble learning model and a wave velocity ellipsoidal physical model, and is equipped with a GPU-accelerated computing unit.

[0065] The adaptive fusion module, with its built-in dynamic weight allocation algorithm and anomaly detection unit, enables optimized fusion of the analytical results from multiple models.

[0066] The 3D reconstruction module, based on octree spatial indexing and B-spline interpolation engine, generates a 3D stress field visualization model with confidence assessment.

[0067] The closed-loop control module connects to the rock drilling rig control system and adjusts drilling parameters and support schemes in real time based on stress analysis results.

[0068] This invention provides a method and system for intelligent analysis of tunnel surrounding rock stress based on drilling parameters. It has the following beneficial effects:

[0069] 1. This invention constructs a high-dimensional feature space by fusing multi-source heterogeneous data of drilling parameters, triaxial wave velocity, and rock mass quality, overcoming the limitations of traditional single-model analysis. The innovative design of the dynamic local coordinate system and mechanical energy calculation unit enables precise quantification of energy characteristics during drilling, providing reliable input for geostress classification.

[0070] 2. This invention leverages the complementary strengths of a physics-driven wave velocity ellipsoid model and a data-driven Stacking ensemble learning model, combined with a dynamic weight allocation strategy, to effectively balance the model's interpretability and generalization ability. Especially in regions with abrupt changes in rock mass integrity, the collaborative analysis of the two models significantly improves the robustness of stress direction identification.

[0071] 3. This invention utilizes octree spatial indexing and B-spline interpolation to model stress fields, achieving a continuous spatial representation of geomechanical parameters. Visualized model overlays and confidence assessments intuitively reveal potential risk areas, providing a scientific basis for real-time adjustment of support parameters.

[0072] 4. This invention forms an intelligent closed loop of "perception-analysis-execution" by directly linking the results of ground stress analysis with the rock drilling rig control system. The system can adaptively adjust parameters such as drilling speed and impact frequency according to the real-time stress state, reducing manual intervention and improving construction safety and efficiency.

[0073] 5. This invention utilizes a difference-driven local recalibration mechanism and a historical data similarity optimization algorithm to achieve rapid correction of model parameters without interrupting operations. This mechanism effectively suppresses error accumulation caused by sensor drift and geological anomalies, ensuring the long-term operational stability of the system. Attached Figure Description

[0074] Figure 1 This is a flowchart of the method of the present invention;

[0075] Figure 2 This is a system architecture diagram of the present invention. Detailed Implementation

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

[0077] Please see Figure 1 This invention provides a method for intelligent analysis of tunnel surrounding rock stress based on drilling parameters. The method includes the following steps:

[0078] S1. Real-time acquisition of drilling parameters, triaxial wave velocity data and rock mass quality parameters through sensor array;

[0079] In this embodiment, a multi-source sensor collaborative acquisition system is used to achieve high-precision synchronous acquisition of drilling parameters, triaxial wave velocity data, and rock mass quality parameters. The sensor group consists of a drilling parameter sensor array mounted on an intelligent rock drilling rig, a triaxial ultrasonic probe group, and a ground-penetrating radar. Each sensor is connected to the data acquisition unit via a CAN bus, and the sampling frequency is not less than 10Hz to ensure dynamic response performance.

[0080] Drilling parameter acquisition: Impact pressure P impReal-time monitoring is achieved using an embedded piezoelectric sensor (range 0-40MPa) in the hydraulic line. The sensor is installed at the oil inlet of the impact piston cylinder, directly reflecting instantaneous pressure fluctuations in the hydraulic system. The rotational torque τ is measured using a flange-type strain torque sensor (accuracy ±0.5%FS), integrated into the output shaft flange of the rotary motor. Data is transmitted wirelessly to eliminate slip ring contact noise. The propulsion speed vector v... drill The measurement is fused between an optical encoder (0.1 mm resolution) and an inertial measurement unit (IMU). The encoder is mounted on the piston rod of the propulsion cylinder, and the IMU provides attitude compensation data to eliminate measurement deviations caused by trolley vibration. Drilling diameter D h With piston rod diameter D r Dynamic calibration was performed using a laser rangefinder (error ±0.2mm). Three sets of rangefinders were arranged at 120° along the circumference of the borehole, and the actual borehole diameter was fitted using the least squares method.

[0081] Triaxial wave velocity data acquisition: The ultrasonic probe group consists of three orthogonally arranged piezoelectric transducers. The distance between the transmitting and receiving probes is fixed at 200mm. Wave propagation time is measured using the pulse-echo method. x-axis wave velocity V x y-axis wave velocity V y Arranged tangentially along the face of the tunnel, with a z-axis wave velocity V z Parallel to the drilling direction, the probe is mounted on the drill pipe guide frame and automatically adjusts its spatial orientation as drilling progresses. Wave velocity calculation employs a cross-correlation algorithm.

[0082]

[0083] Where d is the probe spacing; Δt i To eliminate environmental noise interference, digital filtering (bandpass 50-200kHz) is used to filter the time difference between the transmitted pulse and the first peak of the received signal.

[0084] Rock mass quality parameter analysis: The ground-penetrating radar employs a multi-band synthetic aperture scanning mode (center frequency 800MHz), emitting electromagnetic waves along the borehole axis in a helical trajectory. The rock quality index (RQD) is extracted through the characteristics of the reflected wave signal.

[0085]

[0086] Among them, l intact The length of a complete rock mass segment within the scanned area is identified using the reflected wave amplitude attenuation rate (threshold - 20 dB) and phase continuity criteria; L scan The total scan length was recorded by the radar motion encoder. Data processing employed wavelet transform for noise reduction and Hough transform for detecting the orientation of rock mass fractures.

[0087] Data synchronization and storage: Multi-source data is timestamped using a GPS-synchronized clock (accuracy ±1μs), and the raw data is temporarily stored using a circular buffer queue (capacity 1GB). The data fusion module performs the following operations:

[0088] Impact pressure P imp The torque τ is subjected to a moving average filter (window length 50ms);

[0089] The propulsion velocity vector v drill Transform to the global coordinate system to compensate for changes in the trolley's pose;

[0090] Wave speed data V x V y V z By matching the spatial location of the ground-penetrating radar scan, a mapping relationship is established under the borehole axis coordinate system.

[0091] Preferably, the buffer queue adopts a double buffering mechanism, which performs feature extraction and outlier detection (based on the 3σ criterion) simultaneously during the data writing process to ensure data consistency of subsequent processing modules.

[0092] S2. Establish a dynamic local coordinate system based on the propulsion speed vector direction in the drilling parameters, and calculate the mechanical specific energy according to the impact pressure, rotation torque and rock breaking volume in the drilling parameters.

[0093] In this embodiment, the precise calibration of the drilling spatial reference and the quantitative characterization of rock-breaking energy density are achieved through the construction of a dynamic local coordinate system and the extraction of energy features. The establishment of the dynamic local coordinate system is based on the real-time drilling trajectory and combined with the geometric features of the tunnel face to achieve adaptive adjustment of spatial orientation, providing a unified spatial reference system for subsequent geostress analysis.

[0094] Construction of dynamic local coordinate system: starting from the drilling operation start point Using the spatial reference origin, sub-centimeter-level coordinate calibration is achieved through a trolley positioning system (fusing GNSS and laser SLAM data). The positive z-axis is determined by the propulsion velocity vector v. drill After smoothing by Kalman filtering and normalization, the result is determined as follows:

[0095]

[0096] Point cloud data of the face is acquired via laser scanner along the y-axis. Principal component analysis (PCA) is used to extract the surface normal vector. And after orthogonalization, we obtain:

[0097]

[0098] The x-axis direction is generated from y×z according to the right-hand rule, ultimately forming an orthogonally normalized local coordinate system {x,y,z}. Preferably, the coordinate system update frequency is synchronized with the drilling speed, and coordinate system recalibration is triggered when the change in the thrust direction exceeds 2°.

[0099] Mechanical specific energy calculation: Based on the impact-rotation composite rock breaking mechanism, the calculation of mechanical specific energy (SE) comprehensively represents the energy consumed in the breaking of a unit volume of rock mass. Its expression is as follows:

[0100]

[0101] in, Characterized by the rock-breaking volume in a single drilling cycle, expressed through the borehole diameter D h With drilling depth L pen (Calculated using the cylinder displacement sensor); L is the hydraulic impact piston stroke; P imp Peak impact pressure; D r τ is the diameter of the impact piston rod; τ is the output torque of the rotary motor; v pen For real-time drilling speed; impact rock breaking energy term Reflects the working efficiency of the hydraulic impact system; rotary cutting energy term Characterizing the rotational shear effect of the cutting tool in rock breaking;

[0102] in the formula The term reflects the effect of borehole cross-sectional area on energy distribution, and the denominator v pen Introducing the reciprocal of velocity enhances the energy density sensitivity under low-speed conditions.

[0103] Impact energy term through piston area difference Precisely quantify the equivalent area of ​​the shock wave transmitted to the bottom of the hole, combined with the peak pressure P. imp The impact energy is dynamically calculated based on the stroke L. The rotational energy term uses torque τ and drilling speed v. pen The ratio form effectively reflects the energy deposition of rotary shear power per unit advance. Rock breaking volume V rock As a normalization factor, the total energy is converted into energy consumption per unit volume, giving the mechanical specific energy SE a definite physical dimension (J / m³). 3 This provides standardized input features for subsequent geostress classification.

[0104] Preferably, the mechanical specific energy calculation incorporates a time integral correction during the drilling process:

[0105]

[0106] The integral period T is synchronized with the impact frequency to eliminate calculation errors caused by single impact fluctuations.

[0107] S3. Input the mechanical specific energy and the rock mass quality parameters into the Stacking ensemble learning model to classify the geostress state and output the rock strength stress ratio classification result.

[0108] In this embodiment, a Stacking ensemble learning model is used to achieve multi-classification analysis of geostress state, combining the advantages of physical features and statistical learning to improve the robustness of rock strength-stress ratio classification. The model adopts a two-level heterogeneous architecture, effectively capturing the nonlinear mapping relationship between mechanical specific energy and rock mass quality parameters through the extraction of diverse features from the base model and the linear fusion of the meta-model.

[0109] Stacking ensemble model construction: The first-level base model group includes Support Vector Machine (SVM), K Nearest Neighbors (KNN), and Random Forest (RF), and each model adopts a heterogeneous feature processing strategy:

[0110] SVM model: It uses radial basis function (RBF) to map features to a high-dimensional space. The kernel function parameters are optimized through grid search. When constructing the decision hyperplane, class weights are introduced to balance the sample imbalance problem.

[0111] KNN model: Based on Mahalanobis distance to measure sample similarity, the nearest neighbor number k value is dynamically adjusted to adapt to different rock mass characteristics. The contribution of mechanical specific energy SE and rock quality index RQD is weighted in the distance calculation.

[0112] RF model: Construct multiple decision trees (preferably with a tree depth of no more than 15 layers), set the feature subset sampling ratio to 0.7 to enhance generalization ability, and adopt a hybrid strategy of Gini index and information gain for the splitting criterion.

[0113] The second-level metamodel employs a linear regression model with L2 regularization, whose inputs are the probability vector output from the base models and the three-axis wave velocity interaction term. The feature vector is constructed as follows:

[0114]

[0115] in, This represents the prediction probability of each base model for the i-th type of geostress state; n is the total number of preset geostress categories.

[0116] Rock strength stress ratio calculation:

[0117] The classification results are quantified using the following formula:

[0118]

[0119] Where n is the total number of geostress state categories; p iThe probability of the i-th type of geostress state determined by the model is represented by λ, which is normalized using the Softmax function. i The characteristic coefficients for the corresponding categories are derived from the uniaxial compressive strength σ of the rock mass in historical borehole data. c Compared with the measured value of geostress σ h The regression analysis calibration satisfies

[0120] The input features include normalized mechanical specific energy. Rock quality indicators and the interaction terms of the three-axis wave velocity parameters (Reflects the anisotropic coupling effect of wave velocity).

[0121] Preferably, a feature importance evaluation mechanism is introduced during the model training phase:

[0122] The contribution of SE and RQD is calculated using the permutation feature method, and the input weights are dynamically adjusted.

[0123] The three-axis wave velocity interaction term is transformed using logarithmic transformation. To reduce the impact of dimensional differences (∈=1e-5);

[0124] Category feature coefficient λ i Based on the dynamic updates of regional geological data, coefficient recalibration is triggered when the number of newly added borehole data exceeds the threshold.

[0125] The model is deployed using a hybrid offline-online training mode: in the offline phase, the base model and meta-model parameters are pre-trained using historical data; in the online phase, sliding window incremental learning (window length 50 drill-through cycles) is used, and elastic weight consolidation (EWC) algorithm is used to prevent catastrophic forgetting.

[0126] S4. Based on the geostress state classification results and the triaxial wave velocity data, analyze the three-dimensional wave velocity distribution and construct a wave velocity ellipsoid physical model by fitting the least squares method.

[0127] In this embodiment, a spatial analysis of the three-dimensional wave velocity distribution is achieved by establishing a wave velocity ellipsoid physical model. Combined with geostress classification results and rock mass anisotropy characteristics, the wave velocity response law in the direction of maximum principal stress is accurately characterized. The method is based on the nonlinear coupling relationship between mechanical specific energy and wave velocity, constructing a highly interpretable physical driving model.

[0128] Three-dimensional wave velocity analytical model: for wave velocities V in each orthogonal direction i Establish a quadratic regression relationship between the mechanical specific energy SE and (i∈{x,y,z}):

[0129] V i =α i SE 2 +βi SE+γ i ;

[0130] Where, α i The curvature characteristics reflecting wave velocity changes with energy density are positively correlated with the rock mass's elastic modulus; β i Characterized by the linear response coefficient, affected by Poisson's ratio; γ i The baseline wave velocity term is negatively correlated with the degree of fracture development.

[0131] coefficient α i β i γ i The solution is optimized using the Levenberg-Marquardt algorithm, with initial values ​​set based on a lithological database (e.g., α for granite). x =0.003, β x =0.12, γ x =3200m / s). Preferably, when the local stress is classified as a high stress state, for α i Apply constraint α i ≥0 to ensure the rationality of the physical meaning.

[0132] Wave velocity ellipsoid model construction: ellipsoid parameters are optimized based on the least squares criterion, and the objective function is defined as:

[0133]

[0134] Where A, B, and C represent the semi-axis lengths of the ellipsoid in the x, y, and z axes, respectively, and their physical meanings correspond to the wave velocity distribution characteristics in the direction of maximum principal stress.

[0135] Regularization term λ(A) 2 +B 2 +C 2 To control model complexity and prevent overfitting (preferably, λ = 0.1, determined through cross-validation);

[0136] The data items are presented using squared error to enhance robustness against outliers.

[0137] The ellipsoid parameters are solved using the singular value decomposition (SVD) method:

[0138] Constructing a design matrix With observation vector b = 1;

[0139] Decomposition of M = UΣV using SVD T Find the pseudo-inverse matrix;

[0140] Iteratively update the semi-axis lengths A, B, and C until convergence (residual change rate < 1e). -5 ).

[0141] Preferably, spatial consistency constraints are introduced:

[0142]

[0143] Ensure physical compatibility between the wave velocity gradient field and the geostress classification results. Perform mutual information verification between the principal direction output by the ellipsoidal model and the prediction results of the Stacking model. When the direction deviation exceeds the threshold, trigger local refitting (using the RANSAC algorithm to remove outliers).

[0144] The geostress classification results λ rock Introducing the ellipsoid model as prior knowledge:

[0145] When λ rock >λ threshold When the z-axis semi-axis length C is subjected to the inequality constraint C≥max(A,B), it reflects the consistency between the direction of the maximum principal stress and the drilling direction.

[0146] Interactive items As a coupling factor, it participates in the optimization of ellipsoidal parameters, enhancing the ability to characterize lateral anisotropy.

[0147] The proposed method achieves high-precision analysis of the spatial distribution of wave velocity field through the synergistic optimization of physical constraints and data-driven approaches, providing reliable input features for subsequent multi-model fusion.

[0148] S5. By integrating the analytical results of the wave velocity ellipsoid physical model with the analytical results of the Stacking ensemble learning model, the principal stress direction is calculated using a dynamic weight allocation strategy based on rock mass integrity.

[0149] In this embodiment, a dynamic weight allocation strategy is used to optimize and integrate the analysis results of the physical model and the AI ​​model. By combining the rock mass integrity state and stress field gradient characteristics, the calculation weight of the principal stress direction is adaptively adjusted to improve the analysis accuracy under complex geological conditions.

[0150] The weighting coefficient w of the wave velocity ellipsoid physical model p Determined by a two-factor regulation function:

[0151]

[0152] Among them, w p η is the weighting coefficient of the wave velocity ellipsoid physical model; RQD0 is the preset rock mass integrity benchmark threshold; k is the integrity adjustment factor, which controls the sensitivity of the weights to changes in RQD; η is the stress gradient adjustment factor, which controls the strength of the weights' response to stress abrupt changes. The stress field gradient modulus is calculated using the partial derivatives of the three-dimensional stress tensor; rock mass integrity factor. When the rock quality index RQD is higher than the baseline threshold RQD0, the Sigmoid function output approaches 1, enhancing the weights of the physical model; stress gradient factor Stress field gradient magnitude Calculated using the central difference method:

[0153]

[0154] The adjustment factor η controls the strength of the weights' response to stress abrupt changes (preferably, η = 0.1), and the weights of the AI ​​model are enhanced to cope with nonlinear effects as the gradient increases.

[0155] The final resolved value is calculated through weighted fusion:

[0156] σ final =w p σ phys +(1-w p )σ AI ;

[0157] Where, σ phys The principal stress directions output by the wave velocity ellipsoid physical model are extracted using the cosine of the principal axis directions of the ellipsoid, specifically the direction vector v corresponding to the maximum semi-axis. max =[A,0,0] is transformed to the global coordinate system using a coordinate system rotation matrix; σ AI The principal stress direction angles predicted by the Stacking ensemble learning model are obtained by calculating the dot product of the Softmax probability vector and a preset direction template (such as horizontal structural stress, vertical self-weight stress, etc.).

[0158] Preferably, temporal smoothing is introduced to prevent abrupt changes in weights:

[0159]

[0160] The smoothing coefficient α = 0.2 (corresponding to a time constant of approximately 5 sampling periods). When the ground-penetrating radar detects a sudden change in rock mass integrity (|RQD)... t -RQD t-1 When the weight update rate exceeds 10%, the weight update is temporarily frozen for 3 cycles to ensure data stability.

[0161] The consistency of the model is evaluated by the angle θ between the principal stress direction vectors.

[0162]

[0163] When θ > 15°, confidence assessment is triggered:

[0164] If w p If the value is greater than 0.7, the physical model results will be used first, and the AI ​​model will be marked as anomaly.

[0165] If w p If the value is less than 0.3, initiate the AI ​​model retraining process (using online incremental learning).

[0166] The intermediate results from both models are retained for subsequent difference analysis.

[0167] The method achieves highly reliable analysis of principal stress directions by complementing the advantages of physical interpretability and data-driven approaches, providing accurate input for three-dimensional stress field reconstruction.

[0168] S6. When the relative difference between the analytical results of the wave velocity ellipsoid physical model and the analytical results of the Stacking ensemble learning model exceeds a set threshold, the local recalibration mechanism is triggered and the three-dimensional stress field model is updated.

[0169] In this embodiment, the dynamic updating of the three-dimensional stress field model is achieved through a difference-driven local recalibration mechanism. Combined with octree spatial partitioning and B-spline interpolation algorithm, the adaptability and continuity of the model under complex geological conditions are ensured.

[0170] The relative difference δ is evaluated using normalized error:

[0171]

[0172] When δ > ε (ε = 0.15), it is judged as an abnormal state, triggering a local recalibration mechanism. Preferably, the difference threshold ε is dynamically adjusted according to the geostress classification results: high stress state (λ rock >5) Reduce ε to 0.1 to enhance sensitivity.

[0173] Octree Anomaly Subspace Partitioning: Dynamically generating the anomaly subspace O based on the spatial topological relationship between lidar point cloud data P and real-time drilling trajectory T. s Its boundary conditions are:

[0174]

[0175] in, Boundary values ​​are determined using a point cloud density clustering algorithm.

[0176] Boundary values ​​were determined using the DBSCAN density clustering algorithm:

[0177] Extracting points from the point cloud whose distance from the drilling trajectory T is less than a threshold d th The point set P with a value of 0.5m t ;

[0178] For P t Perform density clustering to identify the centroid c of the largest connected component;

[0179] Centered on c, extend the boundary along the local coordinate system axis by a distance equal to the standard deviation σ of the point cloud distribution.x , σ y , σ z Three times that.

[0180] In subspace O s Internal execution parameter optimization:

[0181] Wave velocity ellipsoid model update: Adjusting coefficient α using constrained particle swarm optimization (PSO) algorithm i ,β i γ i The objective function is:

[0182]

[0183] The constraints include α i ≥0 and the monotonicity of the half-axis length: A≥B≥C.

[0184] Stacking model weight update: Adjust the meta-model weights through online incremental learning, with the loss function being:

[0185]

[0186] Where the reference value σ ref,k Taken from historical similar working conditions (Euclidean distance < 0.2), w is the meta-model weight vector.

[0187] Three-dimensional stress field reconstruction: The stress field is updated using cubic B-spline interpolation (d=3).

[0188] σ(x,y,z)=∑ i,j,k B i,d (x)B j,d (y)B k,d (z)σ ijk ;

[0189] Where σ(x,y,z) is the three-dimensional stress field interpolation function; B i,d The basis function is a d-order B-spline, and the node vector spacing is adaptively adjusted according to the drilling cycle length; σ ijk The stress values ​​at the interpolation control points are given, with an interpolation order d ≥ 3 to ensure stress field continuity. A spatial variogram model is considered.

[0190]

[0191] θ0, θ1, and θ2 are calibrated using maximum likelihood estimation.

[0192] Apply continuity constraints at the subspace boundary:

[0193]

[0194] Where n is the boundary normal direction. Preferably, the virtual control point method is used to insert overlapping nodes at the boundary of adjacent subspaces to ensure the continuity of the first derivative of the stress field.

[0195] The method achieves efficient updating of the three-dimensional stress field through the synergy of local refined modeling and global smooth constraints, providing real-time and reliable geomechanical guidance for tunnel construction.

[0196] Please see Figure 2 The present invention also provides an intelligent analysis system for tunnel surrounding rock stress based on drilling parameters, the system comprising the following modules:

[0197] The multi-source sensing module integrates a drilling-while-drilling sensor array, a three-axis ultrasonic probe, and a ground-penetrating radar to achieve synchronous acquisition of drilling parameters, wave velocity data, and rock mass quality parameters.

[0198] It integrates a high-precision drilling-while-drilling sensor array (including impact pressure, rotational torque, and feed rate vector sensors), a triaxial orthogonal ultrasonic probe group, and a multi-band ground-penetrating radar. The drilling-while-drilling parameter sampling frequency is ≥1kHz. The ultrasonic probes use the pulse-echo method to achieve microsecond-level wave velocity measurement, and the ground-penetrating radar obtains sub-meter-level rock mass quality parameters through synthetic aperture scanning. The data from each sensor are aligned using hardware timestamps and synchronized at the millisecond level using a circular buffer queue.

[0199] The dynamic modeling module includes a local coordinate system generation unit and a mechanical energy calculation unit, which constructs drilling space benchmarks and energy characteristics in real time.

[0200] An integrated local coordinate system dynamic generation algorithm, combined with laser SLAM point cloud registration and principal component analysis, constructs a drilling spatial reference in real time. The mechanical energy specificity calculation unit integrates an impact-rotation energy coupling model, eliminates instantaneous fluctuations through sliding window filtering, and outputs a standardized energy feature vector.

[0201] The hybrid analysis module deploys a Stacking ensemble learning model and a wave velocity ellipsoidal physical model, and is equipped with a GPU-accelerated computing unit.

[0202] A heterogeneous computing architecture is deployed: the CPU runs the nonlinear optimization algorithm for the wave velocity ellipsoid physics model, while the GPU (NVIDIA CUDA acceleration) processes the multi-classification task of the Stacking ensemble learning model in parallel. Model parameters are updated through an online-offline hybrid training mechanism, supporting TensorRT inference optimization.

[0203] The adaptive fusion module has a built-in dynamic weight allocation algorithm and an anomaly detection unit to achieve optimized fusion of the analytical results of multiple models;

[0204] A dual-buffered pipeline design is adopted, with the main pipeline performing dynamic weight fusion calculations and the auxiliary pipeline handling anomaly detection and confidence assessment. An internal state machine enables smooth switching of the fusion strategy, automatically degrading to single-model operation mode when sensor anomalies are detected.

[0205] The 3D reconstruction module, based on octree spatial indexing and B-spline interpolation engine, generates a 3D stress field visualization model with confidence assessment.

[0206] The system utilizes an octree spatial index to achieve efficient management of stress field data, with node resolution dynamically adjusted according to drilling progress (initially 0.5m, minimum 0.1m). The B-spline interpolation engine employs an adaptive node insertion algorithm to automatically densify control points in regions of abrupt stress gradient changes. A visualization model is overlaid with a confidence heatmap, and real-time rendering is achieved using OpenGL.

[0207] The closed-loop control module connects to the rock drilling rig control system and adjusts drilling parameters and support schemes in real time based on stress analysis results.

[0208] Interconnected with the rock drilling rig's PLC system via industrial Ethernet, the control command generation unit maps stress analysis results into drilling parameter adjustment strategies.

[0209] In high-stress areas, the propulsion speed is automatically reduced (≤50% of the rated value) and the impact frequency is increased (≥120% of the reference value);

[0210] Optimize the support scheme for low-integrity areas by dynamically adjusting the anchor spacing and grouting pressure.

[0211] The control cycle is ≤200ms to ensure real-time response to changes in geological conditions.

[0212] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A smart analytical method for tunnel surrounding rock stress based on drilling parameters, characterized in that, The method includes the following steps: S1. Real-time acquisition of drilling parameters, triaxial wave velocity data and rock mass quality parameters through sensor array; S2. Establish a dynamic local coordinate system based on the propulsion speed vector direction in the drilling parameters, and calculate the mechanical specific energy according to the impact pressure, rotation torque and rock breaking volume in the drilling parameters. S3. Input the mechanical specific energy and the rock mass quality parameters into the Stacking ensemble learning model to classify the geostress state and output the rock strength stress ratio classification result. S4. Based on the classification results of the ground stress state and the triaxial wave velocity data, analyze the three-dimensional wave velocity distribution and construct a wave velocity ellipsoid physical model by fitting with the least squares method. S5. By fusing the analytical results of the wave velocity ellipsoid physical model with the analytical results of the Stacking ensemble learning model, the principal stress direction is calculated using a dynamic weight allocation strategy based on rock mass integrity. The dynamic weight allocation strategy adjusts the fusion weights of the physical model and the AI ​​model using both rock mass integrity and stress gradient factors, and its calculation formula is as follows: ； in, These are the weighting coefficients for the wave velocity ellipsoid physical model; These are the rock mass quality parameters; The preset rock mass integrity benchmark threshold; As an integrity adjustment factor, the control weight follows Sensitivity to change; This is a stress gradient adjustment factor that controls the intensity of the weight's response to sudden stress changes; The stress field gradient magnitude is obtained by calculating the partial derivatives of the three-dimensional stress tensor. The final principal stress direction analytical value is calculated by the following formula: ； in, This represents the analytical results of the wave velocity ellipsoid physical model; This represents the predicted output of the Stacking ensemble learning model; S6. When the relative difference between the analytical results of the wave velocity ellipsoid physical model and the analytical results of the Stacking ensemble learning model exceeds a set threshold, a local recalibration mechanism is triggered and the three-dimensional stress field model is updated; wherein, the determination of the relative difference adopts the normalized error evaluation method, the calculation formula of which is: ； in, This represents the analytical results of the wave velocity ellipsoid physical model; This represents the predicted output of the Stacking ensemble learning model; when The local recalibration mechanism is triggered at any time, where The preset difference threshold is used; The recalibration mechanism operates in the abnormal subspace partitioned by the octree. Internally, the parameters of the wave velocity ellipsoid model are updated using a parameter optimization algorithm based on historical data similarity. , , And the weight parameters of the Stacking ensemble learning model.

2. The intelligent analysis method for tunnel surrounding rock stress based on drilling parameters according to claim 1, characterized in that, In step S1: The sensor group includes a drilling parameter sensor, a triaxial ultrasonic probe, and a ground radar installed on the intelligent rock drilling rig. The drilling parameters include impact pressure, which characterizes the operating state of the rock drilling machinery. Rotational torque , propulsion velocity vector Drilling diameter and piston rod structural parameters ; The triaxial wave velocity data is acquired through an ultrasonic array probe arranged around the borehole, including the x-axis wave velocity in an orthogonal coordinate system. y-axis wave velocity z-axis wave velocity , where the z-axis is aligned with the drilling direction; The rock mass quality parameters are rock quality indices obtained through ground-penetrating radar scanning. Its calculation satisfies: ； in, The length of the complete rock mass segment within the scanned area; This represents the total scan length. The drilling parameters, triaxial wave velocity data, and rock mass quality parameters are stored in a buffer queue in a timestamp-aligned manner.

3. The intelligent analysis method for tunnel surrounding rock stress based on drilling parameters according to claim 1, characterized in that, The method for establishing the dynamic local coordinate system in step S2 is as follows: The starting point of the drilling operation is taken as the spatial reference origin. , will propel velocity vector The unit direction vector is defined as the positive z-axis direction, and the geometric surface of the face is obtained by laser scanning. And calculate its normal vector. The x-axis direction is determined by the right-hand rule, which governs the y-axis. Generates in orthogonal directions; The calculation of the mechanical specific energy is based on the combined effect of impact rock-breaking energy and rotary cutting energy, and its expression is: ； in, Characterizing the volume of rock broken in a single drilling cycle, Drilling depth; For hydraulic impact piston stroke; This refers to the peak impact pressure. The borehole diameter; To impact the piston rod diameter; This is the output torque of the rotary motor; This represents the real-time drilling speed.

4. The intelligent analysis method for tunnel surrounding rock stress based on drilling parameters according to claim 3, characterized in that, In the calculation of mechanical specific energy: Impact rock-breaking energy term Reflects the working efficiency of the hydraulic impact system; Rotary cutting energy term Characterizing the rotational shear effect of the cutting tool in rock breaking; Rock breaking volume The calculation of borehole diameter is introduced The square term is used to reflect the effect of the pore wall area on the energy density, the mechanical specific energy. The unit is J / m 3 .

5. The intelligent analysis method for tunnel surrounding rock stress based on drilling parameters according to claim 1, characterized in that, In step S3: The Stacking ensemble learning model adopts a two-level architecture. The first-level base model group includes support vector machine, K-nearest neighbor algorithm and random forest. The second-level meta-model uses linear regression algorithm to weight and fuse the outputs of the base models. The rock strength-stress ratio classification result is calculated using the following formula: ； in, This represents the total number of geostress state categories. The first decision of the model Probability of geomorphic stress state; These are the feature coefficients for the corresponding categories; The characteristic coefficients are pre-calibrated based on the mapping relationship between different rock mass strengths and stress states in historical borehole data, and the input feature vector includes normalized mechanical specific energy. Rock quality indicators and the interaction term of the three-axis wave velocity parameters .

6. The intelligent analysis method for tunnel surrounding rock stress based on drilling parameters according to claim 1, characterized in that, In step S4: The three-dimensional wave velocity analysis is based on the coupling relationship between mechanical specific energy and the anisotropic characteristics of rock mass, and a wave velocity prediction model is established: ； in, , , The correlation coefficient reflects the elastic modulus, Poisson's ratio, and degree of fracture development of the rock mass; Mechanical specific energy; The physical model of the wave velocity ellipsoid is constructed by optimizing the ellipsoid parameters using the least squares method, with the objective function being: ； in, , , These represent the semi-axis lengths of the ellipsoid in the x, y, and z axes, respectively, and their physical meaning corresponds to the wave velocity distribution characteristics in the direction of maximum principal stress. The number of triaxial wave velocity data sampling points participating in the least squares fitting. The sampling point number, , , They represent the first Wave velocity along the x-axis, y-axis, and z-axis of each sampling point in an orthogonal coordinate system; regularization term. Control the complexity of the model to prevent overfitting.

7. The intelligent analysis method for tunnel surrounding rock stress based on drilling parameters according to claim 1, characterized in that, The abnormal subspace The division is dynamically generated based on the spatial topological relationship between lidar point cloud data and real-time drilling trajectory, and its boundary conditions satisfy: ； in, , , , , , Boundary values ​​are determined by point cloud density clustering algorithm; x, y, and z represent the three-dimensional coordinates of point cloud or real-time drilling trajectory points in the dynamic local coordinate system within the anomaly subspace. The three-dimensional stress field model is updated using the B-spline interpolation algorithm, the expression of which is: ； in, This is a three-dimensional stress field interpolation function; , , For the x, y, z directions B-spline basis functions; The stress value at the interpolation control point, and the interpolation order. To ensure the continuity of the stress field, , , For control point indexes.

8. A smart analysis system for tunnel surrounding rock stress based on drilling parameters, applied to the method described in any one of claims 1-7, characterized in that, The system includes the following modules: The multi-source sensing module integrates a drilling-while-drilling sensor array, a three-axis ultrasonic probe, and a ground-penetrating radar to achieve synchronous acquisition of drilling parameters, wave velocity data, and rock mass quality parameters. The dynamic modeling module includes a local coordinate system generation unit and a mechanical energy calculation unit, which constructs drilling space benchmarks and energy characteristics in real time. The hybrid analysis module deploys a Stacking ensemble learning model and a wave velocity ellipsoidal physical model, and is equipped with a GPU-accelerated computing unit. The adaptive fusion module has a built-in dynamic weight allocation algorithm and an anomaly detection unit to achieve optimized fusion of the analytical results of multiple models; The 3D reconstruction module, based on octree spatial indexing and B-spline interpolation engine, generates a 3D stress field visualization model with confidence assessment. The closed-loop control module connects to the rock drilling rig control system and adjusts drilling parameters and support schemes in real time based on stress analysis results.

Images