Methods and systems for monitoring deformation of surrounding rock in coal mine roadways

Abstract

This invention discloses a method and system for monitoring the deformation of surrounding rock in coal mine roadways, relating to the field of geometric measurement technology. The method involves acquiring raw data, removing abrupt changes in displacement data based on principal stress increment data, and calculating pure rheological creep deformation data using the principal stress increment as the integration path. This data is then calculated with the total macroscopic subsidence displacement to obtain the background rigid body settlement data of the strata. Subsequent difference calculations are performed with the axial deformation data to obtain local net compressive deformation data and cooperative compression ratio data. The reflected Bragg wavelength drift data is used to calculate the crack normal width and the dip angle data determined by the principal strain plane orientation of the three-dimensional spatial strain tensor. A spatiotemporal sequence is constructed from the pure rheological creep deformation data, and differential gradient calculations are performed to obtain deformation propagation velocity vector data. Finally, a dimensionless comprehensive deformation characteristic index is obtained. This invention achieves simultaneous measurement of multidimensional geometric deformation of surrounding rock and reconstruction of the spatial distribution of the deformation field.

Description

Technical Field

[0001] This invention relates to the field of geometric measurement technology, specifically to a method and system for monitoring the deformation of surrounding rock in coal mine roadways. Background Technology

[0002] The measurement technology of multidimensional geometric deformation of surrounding rock is an important branch of the field of precision geometric metrology. It aims to accurately reconstruct the spatial distribution of the deformation field of surrounding rock through high-precision measurement of geometric quantities such as displacement, crack width, and dip angle. In recent years, with the advancement of measurement technologies such as fiber optic grating micro-interference and whole-station laser scanning, the measurement of deformation geometry has gradually developed towards multi-source and multi-dimensional spatial information reconstruction.

[0003] Existing measurement technologies mostly measure the spatial linear displacement of the surrounding rock by deploying deep benchmark delamination meters or laser rangefinders to obtain single-dimensional geometric quantities such as the local delamination of the roof. These methods can effectively measure the local delamination dimensions of the rock strata roof under normal working conditions and can output basic measurements such as absolute displacement values ​​or creep rates, possessing a certain geometric measurement capability. However, under extreme working conditions of large deformation accompanied by strong disturbances, the strong impacts generated by violent cutting and blasting introduce instantaneous step jump components into the original displacement measurement signals. Simultaneously, they often cause large-scale and overall subsidence displacement of the measured object. This results in the original dimensional measurements obtained by the rangefinders installed on the roof and infill surface containing both instantaneous elastic jump components and a non-negligible overall subsidence background component, making it difficult to separate and measure the microscopic net compressive deformation geometry and pure fluid deformation geometry of the local key bearing body.

[0004] Furthermore, existing measurement methods (such as CN116697915A) have limited measurement dimensions for crack geometric parameters. Most can only obtain a single linear strain measurement or normal cracking amount, making it difficult to synchronously invert and reconstruct multi-dimensional geometric parameters such as crack normal width and crack dip angle from a one-dimensional fiber strain measurement sequence. At the same time, limited by single-point or single-dimensional measurement methods, existing methods also lack the ability to measure the propagation velocity vector of deformation geometry along spatial distribution, and cannot obtain deformation propagation geometry that has both magnitude and direction indication characteristics. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method and system for monitoring the deformation of surrounding rock in coal mine roadways.

[0006] To achieve the above objectives, the technical solution of the present invention is as follows:

[0007] In a first aspect, the present invention discloses a method for monitoring the deformation of surrounding rock in coal mine roadways, comprising the following steps:

[0008] Acquire the original linear displacement time series data, principal stress increment time series data, reflected Bragg wavelength drift time series data, macroscopic total subsidence linear displacement of the roof surface, and axial original linear deformation of the filling body, which are discretely distributed along the roadway axis.

[0009] Based on the instantaneous loading rate obtained by first-order time-domain differentiation of the principal stress increment time series data, the step abrupt component corresponding to the abnormal loading rate in the original linear displacement time series data is identified and removed to obtain the reference elastic displacement time series sequence; then, the rheological constitutive integral solution of the reference elastic displacement time series sequence is performed using the principal stress increment time series data as the integral stress path to obtain pure rheological creep linear deformation data.

[0010] The pure rheological creep linear deformation data is used as a local shallow delamination reference. The calibrated linear geometric coupling relationship is mapped to local delamination size data. Then, spatial geometric difference operation is performed with the macroscopic total subsidence linear displacement to obtain the formation background rigid body settlement data. Then, the difference operation is performed between the axial original linear deformation and the formation background rigid body settlement data to measure the local net compression deformation data of the infill body. Finally, the cooperative compression ratio data is calculated.

[0011] The reflected Bragg wavelength drift time series data is input into a preset strain space coordinate transformation matrix for inverse transformation calculation to obtain crack normal width data and crack dip angle data determined by the principal strain plane orientation of the three-dimensional spatial geometric strain tensor.

[0012] The spatial-temporal sequence data of surrounding rock geometric deformation is constructed by the pure rheological creep linear deformation data and differential gradient operation is performed to convert it into spatial geometric change rate value and temporal rheological rate value. Then, it is substituted into the characteristic linear velocity solution formula based on the advection model to obtain deformation propagation velocity vector data.

[0013] The cooperative compression ratio data, crack normal width data, crack dip angle data, deformation propagation velocity vector data, and pure rheological creep linear deformation data are extracted, normalized and weighted, and dimensionless comprehensive deformation characteristic indexes are obtained by measurement.

[0014] Secondly, this invention discloses a coal mine roadway surrounding rock deformation monitoring system, which uses the above-mentioned coal mine roadway surrounding rock deformation monitoring method, including:

[0015] The data stream acquisition module is used to acquire the original linear displacement time series data, principal stress increment time series data, reflected Bragg wavelength drift time series data, macroscopic total subsidence linear displacement of the roof surface, and axial original linear deformation of the filling body, which are discretely distributed along the tunnel axis.

[0016] The stripping and solving module is used to identify and remove step abrupt components corresponding to abnormal loading rate moments in the original linear displacement time series data based on the instantaneous loading rate obtained by first-order time-domain differentiation of the principal stress increment time series data, thereby obtaining a reference elastic displacement time series sequence; then, using the principal stress increment time series data as the integral stress path, the reference elastic displacement time series sequence is subjected to rheological constitutive integral solving to obtain pure rheological creep linear deformation data.

[0017] The calculation module is used to use the pure rheological creep linear deformation data as a local shallow delamination reference, map it to local delamination size data through the calibrated linear geometric coupling relationship, and then perform spatial geometric difference operation with the macroscopic total subsidence linear displacement to obtain the formation background rigid body settlement data. Then, the original axial linear deformation variable is performed with the formation background rigid body settlement data to measure the local net compression deformation data of the infill body, and then the cooperative compression ratio data is calculated.

[0018] The crack geometry inversion module is used to input the reflected Bragg wavelength drift time series data into a preset strain space coordinate transformation matrix for inverse transformation calculation to obtain crack normal width data and crack dip angle data determined by the principal strain plane orientation of the three-dimensional spatial geometric strain tensor.

[0019] The solution module is used to construct the spatiotemporal sequence data of surrounding rock geometric deformation through the pure rheological creep linear deformation data and perform differential gradient calculation to convert it into spatial geometric change rate value and temporal rheological rate value, and substitute it into the characteristic linear velocity solution formula based on the advection model to obtain deformation propagation velocity vector data.

[0020] The index evaluation module is used to extract the cooperative compression ratio data, crack normal width data, crack dip angle data, deformation propagation velocity vector data, and pure rheological creep linear deformation data, perform normalized weighted mapping, and measure to obtain a dimensionless comprehensive deformation characteristic index.

[0021] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0022] 1. This invention avoids the erroneous deletion or retention of real creep geometry information by conventional filtering methods, so that the final output of pure rheological creep linear deformation geometry can accurately reflect the long-term stability trend of the surrounding rock under constant load, thus improving the measurement accuracy of the surrounding rock rheological deformation geometry.

[0023] 2. This invention uses the inverse transformation of distributed fiber optic grating wavelength drift data to synchronously invert and reconstruct multidimensional crack geometric parameters such as crack normal width and crack dip angle from a one-dimensional linear strain measurement sequence. At the same time, based on the spatiotemporal sequence differential gradient calculation of pure rheological creep deformation geometry, the measurement dimension and completeness of the geometric characteristics of the surrounding rock deformation field are expanded.

[0024] 3. This invention integrates velocity vector with co-compression ratio, crack geometric parameters, and pure rheological creep deformation data, enabling the system to map multiple monitoring parameters with different physical units and mechanical meanings into a single dimensionless comprehensive deformation characteristic index. This overcomes the deficiency in existing technologies where multiple source geometric quantities lack a unified fusion mechanism, providing a single and clear quantitative basis for roadway support decisions. Attached Figure Description

[0025] The disclosure of this invention is illustrated with reference to the accompanying drawings. It should be understood that the drawings are for illustrative purposes only and are not intended to limit the scope of protection of this invention. In the drawings, the same reference numerals are used to refer to the same parts. Wherein:

[0026] Figure 1 This is a flowchart of the steps of the present invention;

[0027] Figure 2 This is a schematic diagram illustrating the working principle of the present invention;

[0028] Figure 3 This is a system module connection diagram of the present invention;

[0029] Figure 4 This is a flowchart of the system modules of the present invention. Detailed Implementation

[0030] It is readily understood that, based on the technical solution of this invention, those skilled in the art can propose various interchangeable structural methods and implementations without altering the essential spirit of the invention. Therefore, the following detailed embodiments and accompanying drawings are merely illustrative examples of the technical solution of this invention and should not be considered as the entirety of the invention or as limitations or restrictions on the technical solution of this invention.

[0031] In existing techniques for measuring multidimensional geometric deformations of surrounding rock, single-type sensors are often used to measure roof delamination, tunnel convergence, or backfill compression separately, lacking a unified physical correlation model between the measurement data of each geometric quantity. Existing methods struggle to separate the instantaneous elastic jump component triggered by stress fluctuations from the original displacement measurement signal, and the continuous creep deformation component driven by the rheological properties of the rock mass. This results in the measured displacement geometry containing a large amount of noise unrelated to long-term stability. Furthermore, the total roof subsidence geometry is mixed with background rigid settlement and shallow delamination deformation, and the axial compression geometry of the backfill is also coupled with interference from background settlement, making it impossible to independently measure the local net compression deformation geometry of the backfill itself. Furthermore, existing measurement methods have limited measurement dimensions for crack geometry parameters. Most methods can only obtain a single normal crack size, making it difficult to synchronously invert and reconstruct multi-dimensional geometric parameters such as crack normal width and crack dip angle from a one-dimensional fiber strain measurement sequence. They also lack the ability to measure the propagation velocity vector of deformation geometry along the spatial distribution, and cannot obtain deformation propagation geometry that has both magnitude and direction indication characteristics.

[0032] To address the aforementioned issues, this invention eliminates abrupt changes in displacement data based on the instantaneous loading rate of principal stress increment time-series data. It then performs rheological constitutive integral calculations using the principal stress increment as the integration path, thereby separating the pure rheological creep deformation component reflecting the long-term stability of the rock mass from the original displacement signal. Furthermore, through spatial geometric difference operations, the total macroscopic subsidence displacement of the roof is decoupled into the background rigid body settlement and shallow delamination deformation. The background settlement is subtracted from the axial deformation of the infill body to measure the local net compression deformation of the infill body itself, thus calculating the compression ratio index characterizing the coordinated working state of the roof and infill body. Next, the normal fracture width and dip angle are obtained through inverse transformation of reflected Bragg wavelength drift data, and the deformation propagation velocity vector is obtained based on the spatiotemporal sequence differential gradient calculation of the pure rheological creep deformation data. Finally, the multi-source heterogeneous mechanical parameters are normalized and weighted, mapped into a dimensionless comprehensive deformation characteristic index, forming a closed loop for the overall instability risk quantitative assessment of the surrounding rock-infill body coordinated bearing system.

[0033] After introducing the basic concept of the present invention, the embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0034] Example 1:

[0035] To address the problems in existing technologies, such as the obscuring of local pure flow deformation by large overall formation deformation under high-intensity mining scenarios, the difficulty in decoupling and quantifying the coordinated bearing state of the roof and backfill, the inability to synchronously invert multi-dimensional geometric discontinuous displacement parameters of surrounding rock fractures, and the lack of comprehensive instability discrimination based on multi-source heterogeneous deformation information, etc. Figure 1 and Figure 2 As shown, this invention proposes a method for monitoring the deformation of surrounding rock in coal mine roadways, comprising the following steps:

[0036] Step 1: Acquisition and Spatiotemporal Alignment of Multi-Source Heterogeneous Deformation Data Streams:

[0037] The system deploys multi-type sensor arrays along the axial direction of the coal mine roadway, forming a multi-physics field raw data acquisition layer. Several discrete displacement sensing nodes are deployed along the roadway roof and sides at preset intervals Δs (typically 0.5m to 2m). These nodes preferably employ laser ranging units or fiber optic displacement sensors, outputting raw linear displacement time-series data along the roadway axial direction. Where s represents the discrete spatial coordinates along the tunnel axis and t represents the sampling time; vibrating wire or fiber optic grating stress sensors are deployed at key stress bearing points in the surrounding rock to output principal stress increment time-series data Δσ(t), with a sampling frequency preferably not less than 1Hz; a distributed fiber optic grating (FBG) bundle array is laid in a pre-drilled hydraulic fracturing borehole, the fiber optic grating bundles penetrating the measured fracture area in a spatially orthogonal or cross-oriented manner to output reflected Bragg wavelength drift time-series data Δλ(t); a total-station laser scanner or deep base point delamination meter is deployed on the roof surface to output the macroscopic total subsidence linear displacement of the roof surface. Displacement measuring devices are installed on the surface of the filling material to output the original axial linear deformation of the filling material. .

[0038] All sensor nodes are uniformly connected to the data acquisition front end, and microsecond-level timestamp alignment is achieved through the IEEE 1588 Precision Clock Synchronization Protocol (PTP). The spatial position of each sensor node is based on the axial mileage of the tunnel as a unified spatial reference coordinate (the origin is taken as a fixed measuring station, and the direction of excavation is positive). During installation, the mileage coordinates of each displacement node, stress sensor, separation meter, filling displacement gauge, and fiber optic grating area are measured using a total station and registered in the database. Subsequent spatial registration uses the registered mileage of each sensor as an index to uniformly map it to a regular grid based on the sampling position of pure rheological creep deformation data, thereby ensuring that multi-source heterogeneous data are fused in the same spatiotemporal coordinate system.

[0039] Step 2: Elimination of mutation components and solution of rheological constitutive integrals:

[0040] To address the issue that the intense impact disturbances generated by the violent cutting of coal mining machines and the blasting of large-section arches introduce instantaneous step-jump deformations into the original linear displacement time-series data, this invention employs a solution strategy combining disturbance identification based on the instantaneous loading rate of principal stress increments with rheological constitutive integrals. This strategy extracts pure rheological creep linear deformation data that truly characterizes the rheological properties of the surrounding rock from the original linear displacement time-series data. The specific steps are as follows:

[0041] First, extract the principal stress increment time series data. Perform a first-order time-domain derivative operation on it to obtain the instantaneous stress loading rate. : ; Instantaneous stress loading rate Compared with the preset rate of change threshold Numerical comparisons are performed to identify abnormal instantaneous moments that exceed a preset rate of change threshold. And lock onto the abnormal instantaneous time in the original linear displacement time series data u(s,t). Corresponding step jump deformation size Preset rate of change threshold The value is determined based on the physical and mechanical parameters of the surrounding rock on site, with a typical range of 0.5 MPa / s to 2.0 MPa / s.

[0042] Identified step jump deformation size The components were defined as abrupt changes and removed; the resulting displacement sequence was denoted as the reference elastic displacement time series. :

[0043] ;

[0044] H(·) is the Heaviside step function, used to maintain baseline continuity after the removal of points.

[0045] Subsequently, using principal stress increment time series data As a stress-rheological integration path, the reference elastic displacement time series Perform time-domain hardening integral operations to generate pure rheological creep linear deformation data using Stieltjes convolution. That is, for the infinitesimal element of stress increment Perform convolution integral:

[0046] ;

[0047] in, The pre-calibrated constitutive compliance kernel function of the surrounding rock reflects the rheological response characteristics of the surrounding rock under different stress levels and loading durations, using time difference... The independent variable represents the memory effect of the surrounding rock rheological response on the stress loading history. The compliance kernel function J parameter was obtained by laboratory uniaxial creep tests.

[0048] The kernel function of the surrounding rock rheological constitutive compliance is constructed based on the Burgers rheological constitutive model. The Burgers model consists of a Maxwell element and a Kelvin-Voigt element connected in series, used to simultaneously characterize the three types of responses of the surrounding rock: instantaneous elastic, steady-state viscous flow, and recoverable viscoelastic creep. Under unit step stress, the explicit expression of its creep compliance function J(t) is:

[0049] ;

[0050] In the formula, This is the instantaneous elastic modulus of the Maxwell element (characterizing the elastic strain response at the instant of loading, typically ranging from 1 to 20 GPa). The viscosity coefficient of the Maxwell element (characterizing irreversible viscous flow during the steady-state creep stage, typically ranging from 1 × 10⁻⁶). 13 ~1×10 15 Pa·s); The viscoelastic modulus of the Kelvin-Voigt element (characterizing the recoverable deformation amplitude during the transition creep stage, typically ranging from 0.5 to 10 GPa). The viscosity coefficient of the Kelvin-Voigt unit (characterizing the time characteristics of transition creep, typically ranging from 1 × 10⁻⁶). 12 ~1×10 14 Pa·s); t is the loading duration; τ is the stress loading time. These four parameters constitute the set of surrounding rock rheological constitutive parameters Θ={ , , , }

[0051] Furthermore, when the stress Δσ(τ) is continuously differentiable along the loading path, according to the Boltzmann principle of linear superposition, The Stieltjes convolution integral is equivalent to the following Riemann convolution integral form:

[0052] ;

[0053] In the formula, the lower limit of integration Take the reference elastic displacement time series after removing abrupt changes. The initial sampling time corresponds to the initial conditions. That is, the monitoring start time is taken as the common zero reference for strain and stress; the upper limit of integration t is taken as the current sampling time.

[0054] When the surrounding rock exhibits characteristics dominated by recoverable viscoelastic creep and negligible steady-state viscous flow, the generalized Kelvin-Voigt model can be preferred for constructing the compliance kernel function, which consists of N Kelvin-Voigt elements connected in series, with the number of elements N preferably being 2 to 3.

[0055] ;

[0056] in , Let θ be the viscoelastic modulus and viscosity coefficient of the i-th element, and let Θ = { , The parameter set Θ is obtained by conducting laboratory uniaxial or triaxial graded creep tests on standard rock samples collected from the surrounding rock at the site. The rock samples are graded and loaded to 50%–80% of the engineering stress level. A constant stress is maintained at each stress level, and the creep curves showing the change in axial strain over time are recorded. Then, using the compliance function J(t) as the fitting model, a nonlinear least squares method (such as the Levenberg-Marquardt algorithm) is employed to fit the creep curves, solving for the optimal values ​​of each parameter in Θ. A typical fitting correlation coefficient R² ≥ 0.95 is achieved, and the parameter calibration error is controlled within ±5%. The calibrated parameter set Θ is pre-stored in a calibration parameter library.

[0057] This application uses time-domain hardening integral calculation to accumulate the internal strain of the surrounding rock based on the current stress path. Simultaneously, it eliminates high-frequency symmetrical mechanical vibration noise (such as zero-mean white noise introduced by equipment vibration and electromagnetic interference) with zero mean by mean filtering of adjacent sampling points, generating pure rheological creep linear deformation data to characterize the true rheological response of the surrounding rock. Pure rheological creep linear deformation data eliminates instantaneous elastic jumps and high-frequency vibration interference, retaining only the continuous creep rheological information determined by the rock mass's own properties, providing physically realistic input for subsequent collaborative compression ratio calculation and spatiotemporal sequence construction.

[0058] Step 3: Calculation of rigid body sedimentation coupling and synergistic compressibility ratio in the stratigraphic background:

[0059] The pure rheological creep linear deformation data output from step two As a local shallow delamination reference, it maps to the local delamination size data of shallow unstable rock layers deviating from the deep stable anchoring point. The physical basis lies in the following: Under the stress boundary conditions of the tunnel roof, creep along the tunnel axis is mainly induced by the bending and subsidence of the roof strata; according to the Poisson effect in elasticity, when vertical delamination of the roof strata occurs, a proportional axial strain will be generated. Therefore, there is a linear proportional relationship between the vertical delamination and the axial creep, determined by both the Poisson ratio of the strata and the geometric constraints of the roof bending, rather than their numerical values ​​being equal. Based on this, this method is applied to the monitored characteristic measurement points. Introduce calibrated linear geometric coupling relationships at this point:

[0060] ;

[0061] In the formula, The geometric coupling calibration coefficient is a comprehensive reflection of the geometric constraints of the Poisson's ratio of the rock strata and the bending and settlement of the roof. The calibration method is as follows: at the feature measurement points Vertical delamination meters and axial displacement sensors are deployed simultaneously to collect synchronous data for at least one stable production cycle. The proportionality coefficient between the vertical delamination amount and the axial creep variable is determined using least squares regression. And cross-validate the data using data from independent time periods, requiring a correlation coefficient R. 2 Not less than 0.9; when independent vertical delamination calibration data are lacking on site, take... =1 is used as a conservative engineering approximation (at this time) The applicability of this approximation is periodically verified by on-site vertical separation meter measurements. Recalibration is performed when the deviation exceeds a preset limit. .

[0062] To address the problem of the intertwined overall ground subsidence and local net compression height of the bearing body under high-intensity mining scenarios, which traditional methods cannot decouple, this invention separates the rigid settlement component of the ground background from the total macroscopic subsidence displacement. This allows for precise measurement of the actual pressure change of the backfill material itself, and the linear displacement of the total macroscopic subsidence on the roof surface is then calculated. Subtract local delamination dimensions The background rigid body settlement data of the strata were obtained by difference decoupling. :

[0063] ;

[0064] Stratigraphic background rigid body settlement data The rigid translation component, which characterizes the overall strata due to the movement of the overlying strata in the goaf, is unrelated to the local bearing capacity of the filling body and should be removed as an interference in subsequent calculations.

[0065] Because the axial deformation of the filling body is also coupled with the background rigid body settlement, the original linear axial deformation of the filling body is... Rigid body settlement data with stratigraphic background By performing interpolation calculations, background overall settlement interference is completely eliminated, and local net compression deformation data characterizing the actual change in pressure of the infill material itself are obtained. :

[0066] ;

[0067] Finally, the local net compression deformation data As a numerator, the local delamination size data As the denominator, a division ratio calculation is performed, outputting the cooperative compression ratio data used to quantitatively characterize the reverse support constraint state between the top surface and the infill. :

[0068] ;

[0069] Collaborative compression ratio data For example: when the delamination dimension of the top slab increases, if the filling material can simultaneously provide a proportional compression response, then... Approaching the preset support design collaborative benchmark value (Typical values ​​range from 0.6 to 0.9), indicating a good reverse support constraint between the top plate and the filling material; when A significant deviation from this benchmark value indicates an abnormality in the stiffness matching between the two, suggesting a risk of local instability. Preferably, when When this occurs, the system outputs an early warning message indicating abnormal stiffness of the support system.

[0070] Step 4: Inversion of geometric parameters of cracks in distributed fiber gratings:

[0071] To simultaneously acquire the normal crack size and dip angle geometric parameters of surrounding rock fractures, this invention utilizes the reflected Bragg wavelength drift time-series data output from a distributed fiber optic grating array laid within a hydraulic fracturing borehole. Through a preset strain space coordinate transformation matrix, an inverse transformation is performed to reconstruct the fracture geometry parameters from one-dimensional fiber strain. The specific steps are as follows:

[0072] First, the reflected Bragg wavelength drift time series data Δλ(t) is combined with the preset photoelectric strain sensitivity coefficient. Linear mapping was performed to obtain a one-dimensional multi-orientation line strain measurement sequence along the fiber axis. :

[0073] ;

[0074] in, This refers to the reflected Bragg wavelength drift timing data obtained in step one, which is the offset of the reflected wavelength of the fiber grating relative to the initial center wavelength. This is the center wavelength of the fiber Bragg grating (typical value 1550nm). The photoelectric strain sensitivity coefficient (typically around 0.78, corresponding to approximately 1.2 pm / με at a wavelength of 1550 nm). One-dimensional multi-orientation line strain measurement sequence. The strain vector is formed by measuring multiple fiber grating bundles distributed at different spatial orientations (such as 0°, 45°, 90°, and 135°). .

[0075] Subsequently, the one-dimensional multi-orientation line strain measurement sequence is input into the strain space coordinate transformation matrix. The inverse transformation matrix multiplication operation is performed to invert and reconstruct the three-dimensional spatial geometric strain tensor measurement data. Strain space coordinate transformation matrix Based on the spatially orthogonal or intersecting geometric topology of the distributed fiber grating bundles laid within the hydraulic fracturing borehole, a geometrically compatible equation in solid mechanics is constructed, whose elements are quadratic combinations of the direction cosines between each oriented fiber and the principal strain axis in three-dimensional space. The inverse transformation solution formula is as follows:

[0076] ;

[0077] in, It is the inverse matrix or generalized inverse of the strain space coordinate transformation matrix (the least squares generalized inverse is used when the orientation number is greater than 3).

[0078] Introducing the discontinuous boundary condition of crack displacement from fracture mechanics, and measuring data from three-dimensional geometric strain tensors. The amount of open geometric discontinuity displacement was measured by peeling. Shear-type geometric discontinuity displacement The integral of strain along the crack normal (i.e., the direction of the crack surface normal vector) directly corresponds to the opening crack width, while the integral of shear strain along the crack surface tangentially corresponds to the shear slip. Crack normal width data. The formula is:

[0079] ;

[0080] Crack dip angle data From the three-dimensional spatial geometric strain tensor The normal orientation of the principal strain plane is determined. For Perform eigenvalue decomposition and extract the eigenvectors corresponding to the crack surface normal. (Determined by the direction of minimum principal strain), let the unit axial vector of the tunnel / rock mass be... (Taking the x-axis), the dip angle of the fracture surface relative to the axis of the surrounding rock being measured is:

[0081] ;

[0082] Furthermore, the quantity originally defined as the ratio of opening-type and shear-type discontinuous displacements is renamed the crack displacement mixing angle. As an independent parameter, it characterizes the degree of mixing of failure modes during the transition from opening-type to shear-type cracks:

[0083] ;

[0084] Where L is the characteristic length of the fiber grating bundle spanning the crack. This represents the strain distribution along the normal direction. Characterizes the linear size of the crack in the direction perpendicular to the crack surface (typically within the range of 0.05 mm to 10 mm in engineering). The dip angle of the principal strain plane of the fracture surface relative to the axis of the surrounding rock being measured reflects the shear propagation characteristics of the fracture (the value ranges from 0° to 90°).

[0085] Preferably, before constructing the spatiotemporal sequence data of surrounding rock geometric deformation, a step to enhance the spatial consistency of deformation data is also included:

[0086] The system acquires the total macroscopic linear displacement of the top surface. The original linear deformation of the filling body along its axial direction Spatial coordinate layout, using pure rheological creep linear deformation data Using the reference grid as the sampling location, spatial resolution interpolation resampling (preferably cubic spline interpolation or B-spline interpolation) is performed on the total macroscopic subsidence linear displacement and the original axial linear deformation to achieve grid dimension alignment of the multi-source heterogeneous deformation data, and output the spatially resolution aligned multi-source deformation dataset. .

[0087] Simultaneously, the system collects local temperature increment data ΔT(t) caused by the heat release from the hydration of underground filling (obtained through fiber optic grating temperature sensors or platinum resistance temperature sensors deployed within or adjacent to the filling body), and calculates the data based on a preset fiber thermal expansion correction coefficient. For one-dimensional multi-orientation line strain sequences Perform temperature compensation calculations:

[0088] ;

[0089] in, The equivalent thermal expansion coefficient of the fiber Bragg grating (typically around 6.5 × 10⁻⁶). -6 / ℃). Temperature-compensated one-dimensional multi-orientation line strain sequence Replace the original one-dimensional multi-orientation strain sequence input strain space coordinate transformation matrix An inverse transformation operation is performed to eliminate the influence of the hydration heat temperature effect on the accuracy of crack geometry parameter inversion.

[0090] Step 5: Spatiotemporal sequence construction and velocity vector calculation for deformation front kinematics propagation:

[0091] To achieve advanced localization of the propagation front of the surrounding rock fracture surface within the tunnel space, this invention constructs a spatiotemporal sequence containing two-dimensional geometric features along both spatial and temporal axes based on pure rheological creep linear deformation data. By simultaneously solving differential gradient calculations and characteristic line equations, it outputs deformation propagation velocity vector data with directional indication characteristics. The specific steps are as follows:

[0092] Pure rheological creep linear deformation data generated at multiple discretely distributed measuring points along the tunnel axis. Perform spatial coordinate alignment (where (representing the axial coordinate of the i-th measuring point), and introducing the unified high-precision clock stamp from step one. Construct and generate a spatiotemporal sequence data matrix of surrounding rock geometric deformation containing two-dimensional geometric features of spatial axis s and time axis t. .

[0093] The preset spatiotemporal grid difference operator is invoked as the execution carrier for the difference gradient operation, and the spatiotemporal sequence data matrix of surrounding rock geometric deformation is processed. Perform adjacent discrete grid difference calculations along the spatial axis s of the tunnel, and output the spatial geometric change rate value. :

[0094] ;

[0095] Synchronously perform time-domain adjacent time step differential calculations along the time axis t, and output the time rheometry rate value. :

[0096] ;

[0097] Where Δs is the spatial interval between adjacent measurement points and Δt is the sampling interval between adjacent time points, the central difference scheme is adopted to improve numerical stability.

[0098] Furthermore, this application approximates the pure rheological creep deformation field of the surrounding rock as a kinematic wave transported unidirectionally along the tunnel axis, and its governing equation is the advection equation:

[0099] ;

[0100] Based on the characteristic line theory of mathematical physics equations, and assuming that u remains constant along the characteristic line, the kinematic propagation velocity of the deformation equivalent front along the tunnel axis can be obtained. :

[0101] ;

[0102] The negative sign in the formula indicates that the propagation direction is opposite to the direction of the spatial deformation gradient (deformation spreads from the high-value region to the low-value region). It represents the ratio of the time-varying rate of rheology to the spatial geometric rate of change.

[0103] Deformation propagation velocity vector It represents the kinematic spread speed of the deforming forward, used to anticipate and locate the breaking forward. Direction determination needs to consider both time-varying rate and rheological rate. With spatial geometric change rate A combination of symbols. Specifically:

[0104] when When the velocity is >0 (surrounding rock creep continues to increase, which is the most common precursor to instability in engineering), the velocity direction is from... The sign of the symbol is determined by the negative sign in the formula:

[0105] like >0 (deformation increases along the roadway axis), then <0 indicates the negative axial direction;

[0106] like <0, then >0 indicates the positive axial direction;

[0107] when The arithmetic mean of the effective velocity values ​​in the neighborhood of the current measurement point (within 1 to 2 spatial steps before and after) is taken as the deformation propagation velocity at that point. ε is a preset numerical stability threshold, which is 5 × 10⁻⁶. -4 .

[0108] when When the value is less than 0 (rare conditions, such as short-term elastic rebound of surrounding rock or sensor malfunction), the direction determination rule is reversed. This system defaults to using... >0 provides direction output for the engineering baseline operating condition, and in When the absolute value is less than the preset noise threshold (typical value 0.01 mm / h), the direction sign remains unchanged to enhance the stability of the engineering criterion.

[0109] Modulus of deformation propagation velocity vector data Characterizing the rate of deformation propagation, for example A value >0.5 m / h indicates a relatively rapid deformation propagation trend.

[0110] Step Six: Multi-source parameter normalization weighted mapping and comprehensive deformation characteristic index measurement:

[0111] To unify the mapping of various deformation monitoring parameters with different physical units and mechanical meanings into a single dimensionless index, this invention performs a normalized weighted mapping operation based on the design limit geometric feature matrix and normalized weighted expert vectors, outputting a dimensionless comprehensive deformation feature index that can be directly used for instability detection. The specific steps are as follows:

[0112] First, the system retrieves the preset design limit geometric feature matrix. Design the limit geometric feature matrix Including the limit deformation threshold Limiting compression ratio threshold Limiting propagation speed threshold and the limit seam width threshold This constitutes the four-dimensional design limit vector:

[0113] ;

[0114] The values ​​of each limit threshold are pre-calibrated based on the on-site surrounding rock mechanical parameters, tunnel support design standards, and historical disaster statistics. Typical value ranges are provided for reference. The value range is 50mm to 200mm (adjusted according to the specific tunnel cross-section dimensions). The value range is 1.0 to 1.5. The value range is 0.5 m / h to 2.0 m / h. The value range is 5mm to 20mm.

[0115] Subsequently, the preset tilt angle weighting correction function is obtained. (·), This will display the crack dip angle data. Input the tilt angle weight correction function to calculate the seam width weight correction coefficient. :

[0116] ;

[0117] Where α is the dip angle influence coefficient, with a value ranging from 0.2 to 0.5. The physical meaning of the dip angle weighting correction function is: when the crack dip angle is closer to 45°, the influence of shear-type discontinuous displacement on the instability of the surrounding rock is more significant, and the crack width needs to be weighted and amplified for correction.

[0118] Crack normal width data Multiply by seam width weighting factor Obtain the weighted corrected crack normal width data. :

[0119] ;

[0120] Pure rheological creep linear deformation data (Take the maximum value or representative measurement point value along the roadway axis), cooperative compression ratio data Magnitude in deformation propagation velocity vector data and the weighted corrected crack normal width data , respectively with the corresponding limiting deformation threshold Limiting compression ratio threshold Limiting propagation speed threshold and the limit seam width threshold Perform matrix division normalization to generate multidimensional dimensionless geometrically degenerate eigenvectors. :

[0121] ;

[0122] Finally, extract the pre-defined normalized weighted expert vector. (satisfies Σ) =1, i=1,2,3,4), the multidimensional dimensionless geometrically degenerate eigenvectors Perform a matrix inner product multiplication and addition mapping with the normalized weighted expert vector W to output a dimensionless comprehensive deformation feature index. :

[0123] ;

[0124] Among them, the typical values ​​of each component of the normalized weighted expert vector W are: =0.35 (shape weight) =0.25 (Cooperative compression ratio weight) =0.20 (propagation speed weight) =0.20 (joint width weight), the specific value is determined based on the characteristics of the surrounding rock in the mining area and the experience of support design experts. Dimensionless comprehensive deformation characteristic index It is a single dimensionless value, directly related to the overall stability margin of the surrounding rock-fill bearing system. This indicates that the system is in a safe state. This indicates that the system is approaching or exceeding its design limits, and corresponding control measures need to be initiated.

[0125] Furthermore, a verification mechanism is set for the limit threshold and weight vector: using historical disaster samples and stable operating condition samples as calibration and verification sets, the limit threshold and weight vector are back-checked using the hold-out method. When the false alarm rate / false alarm rate of the back-check exceeds the preset limit, recalibration is triggered; after the weight vector is updated, it must satisfy the constraint that each component is non-negative and the sum is 1.

[0126] Step 7: Issuance and execution of physical closed-loop control commands:

[0127] The dimensionless comprehensive deformation characteristic index was obtained by measurement. Subsequently, the system further executes physical closed-loop control commands, achieving an end-to-end closed loop from data measurement to engineering response. The specific steps are as follows:

[0128] The system will integrate dimensionless deformation characteristic indices (t) and the preset safety critical control threshold (Typical value 1.0) Performs numerical comparison logic operations to identify... The deformation exceeding the limit event.

[0129] In response to deformation exceeding the limit event, the system extracts deformation propagation velocity vector data. Spatial direction symbols included And combined with the velocity vector magnitude Locking the dangerous space boundary pointed to by the deformation propagation front within the alleyway space. Specifically, based on the current location of the measurement point with the greatest deformation. Starting from, along The indicated direction is extrapolated to a preset warning distance. (Typical values ​​range from 5m to 20m), located at arrive The tunnel space section within the range is defined as the hazardous space boundary. .

[0130] Based on the dimensionless comprehensive deformation characteristic index The difference from the safety critical control threshold The system sends a reverse stiffness compensation control command to the inverter drive of the flexible mold pumping grouting filling equipment deployed within the boundary of the hazardous space, and sends a secondary compensation grouting command to control the grouting pump to compensate for the grouting pressure. :

[0131] ;

[0132] in, The initial reference pumping pressure for the filling equipment (determined by the equipment's rated parameters) is β, which is the pressure gain coefficient (typically ranging from 0.3 to 1.0). Reverse stiffness compensation control commands are directly sent to the inverter drive of the filling equipment via an industrial fieldbus (such as PROFIBUS or Modbus-TCP), achieving an end-to-end physical execution closed loop.

[0133] Step 8: Environmental Adaptive Online Calibration

[0134] To address the problem of time-varying drift of surrounding rock physical and mechanical parameters during mining stress release, leading to a gradual mismatch between the initially calibrated constitutive model and normalized weighted expert vectors and actual field conditions, this invention introduces an environmentally adaptive online calibration mechanism after the physical closed-loop control command is issued to dynamically correct key system parameters. The specific steps are as follows:

[0135] After issuing the reverse stiffness compensation control command, the system continuously acquires the boundary of the dangerous space. The original linear displacement time series data u(s,t) and principal stress increment time series data were reacquired within the time frame. , as a closed-loop feedback input.

[0136] Based on the principal stress increment time series data at the current moment Call the pre-calibrated constitutive model of surrounding rock creep. And its parameter set Θ, calculate the actual stress increment rate of change after reverse stiffness compensation. And substitute it into the constitutive model of the surrounding rock creep. Perform time-domain integration to generate the expected rheological displacement. :

[0137] ;

[0138] Subsequently, the measured displacement increment of the original linear displacement time series data u(s,t) at the current time step is calculated. And calculate its relationship with the expected rheological displacement. The actual deformation reconstruction residual r(t) between:

[0139] ;

[0140] Perform a first-order difference operation with respect to the time dimension on the actual deformation reconstruction residual r(t) to obtain the residual evolution rate value. :

[0141] ;

[0142] When the residual evolution rate value Greater than the preset convergence threshold When the value is typically 0.5 mm / h, it is determined that the surrounding rock is experiencing nonlinear accelerated failure disturbance, and the constitutive model deviates significantly from the actual response. In this case, the system retrieves the current principal stress increment time-series data Δσ(t) as a disturbance correction factor. ( Using the reference stress increment value (the statistical mean of the principal stress increments during the first stable production cycle after system startup, such as within 72 hours), the normalized weighted expert vector W is updated using linear adaptive gain compensation.

[0143] ;

[0144] Wherein, γ is the adaptive learning rate (typically ranging from 0.01 to 0.05). This represents the baseline feature vector recorded during the system startup phase under stable operating conditions. The updated normalized weighted expert vector... Ensure the updated weight vector The sum of each component is 1, and it is written back to the normalized weighted mapping operation stage, so that the measurement results of the subsequent comprehensive deformation characteristic index are more in line with the real mechanical response under time-varying working conditions.

[0145] Through the system integration of the above eight steps, this invention achieves physical separability of elastic jumps and rheological responses in the original displacement data by mutation elimination and rheological integral calculation; the rigid sedimentation decoupling of the formation background eliminates the masking of local net deformation by the overall large deformation; the distributed fiber grating crack geometric inversion synchronously outputs opening and shear discontinuous displacement parameters; the multi-source parameter normalized weighted mapping realizes the instability discrimination of heterogeneous multiphysics data; and the physical closed-loop control command issuance and environmental adaptive online calibration further form a closed loop from measurement to response and from response to self-calibration. The technical effects of each module are mutually enhanced, constituting a complete multiphysics surrounding rock stability monitoring system.

[0146] In a preferred embodiment of the present invention, the parameter configuration and engineering implementation details for mutation component elimination and rheological constitutive integral solution are as follows:

[0147] The sampling frequency of the principal stress increment sensor is preferably no less than 10Hz to ensure sufficient time resolution for instantaneous stress jumps caused by coal mining machine cutting and blasting impact. The derivative of the instantaneous stress loading rate σ̇(t) is calculated using a forward difference or central difference scheme, and a low-pass filter with a cutoff frequency of 5Hz is pre-applied to prevent high-frequency quantization noise from being amplified during the derivative calculation. A preset rate of change threshold is also included. The calibration method is as follows: Collect continuous 24-hour principal stress increment time-series data under normal tunnel excavation and production conditions, statistically analyze the instantaneous loading rate distribution, and take the upper limit of its 99.7% confidence interval (i.e., mean + 3 standard deviation) as the baseline. When the on-site working conditions change significantly, the data can be updated on a sliding basis based on the data from the most recent 7 days.

[0148] Mutation component The identification uses a front and back window comparison method: Centered on the forward window With backward window ( (Typical values ​​are taken as 0.5s to 2s), and the average displacement within the two windows is calculated respectively. and ,but .

[0149] In another preferred embodiment of the present invention, the engineering parameters and anomaly determination logic for the coupling of formation background rigid body sedimentation and the calculation of synergistic compressibility ratio are as follows:

[0150] Macroscopic total subsidence linear displacement Measurements are preferably performed using a total-station laser scanner or a deep benchmark delamination instrument, with measurement points located at least 30m away from the working face; the original linear deformation of the axial direction of the filling body. The measurement is preferably performed using a fiber optic strain gauge or a micro-displacement sensor, with a measurement accuracy better than 0.1 mm.

[0151] Collaborative compression ratio data The engineering anomaly determination adopts a three-level classification mechanism: when In When, it is judged as a normal state; when In When this occurs, it is determined to be a warning state, triggering on-site inspection; when or If the situation is deemed dangerous, an audible and visual alarm will be triggered immediately, and the mine pressure management and dispatch center will be notified.

[0152] In another preferred embodiment of the present invention, the hardware configuration for distributed fiber Bragg grating gap geometry inversion is as follows:

[0153] The preferred borehole angle for hydraulic fracturing is perpendicular to the expected main fracture surface, with a depth range of 5m to 15m and a diameter of Φ32mm to Φ50mm. The distributed fiber grating array consists of four independent optical fibers, arranged along four spatial orientations: 0°, 45°, 90°, and 135°. The spacing between adjacent grating regions is preferably 10mm to 50mm, and the total number of grating regions in a single borehole is no less than 100. The fiber demodulator preferably employs wavelength scanning demodulation technology, with a wavelength measurement accuracy better than ±1pm and a corresponding strain measurement accuracy better than ±1με.

[0154] For a 4-orientation fiber bundle, the strain space coordinate transformation matrix is... It is a 4×3 dimensional matrix, where each row of elements is determined by the corresponding fiber orientation angle. With the three-dimensional principal strain axis The quadratic combination of the direction cosines between them forms the following specific form:

[0155] ;

[0156] matrix The inverse operation employs the Moore-Penrose generalized inverse to address the least-squares solution for overdetermined systems. Preferably, when the residual norm exceeds a preset threshold (typically 5με), it is determined that there is an anomaly in the fiber grating bundle layout or that the crack geometry model deviates from the single-plane assumption, requiring the initiation of a multi-plane segmented inversion mechanism: the measurement interval is divided into several local segments along the fiber axis, and a local strain space coordinate transformation matrix is ​​constructed for each segment. The inversion is performed independently, and then fused by geometric constraints of adjacent segments to improve the inversion accuracy under complex crack morphology.

[0157] Fiber thermal expansion correction coefficient The calibration is usually obtained by temperature scanning calibration of the unstressed free segment of the fiber grating in a temperature-controlled chamber. The calibration temperature range covers 0℃ to 60℃, and the calibration residual is controlled within ±0.5με / ℃.

[0158] In another preferred embodiment of the present invention, the engineering parameters for constructing the spatiotemporal sequence and calculating the deformation propagation velocity vector are as follows:

[0159] Spatiotemporal sequence data matrix of surrounding rock geometric deformation The optimal spatial grid density is as follows: the spatial interval Δs between adjacent measuring points ranges from 0.5m to 2m, the time sampling interval Δt ranges from 1s to 10s, and the spatiotemporal grid resolution should match the physical spatiotemporal scale of the deformation propagation of interest. The preferred spatiotemporal grid difference operator is a fourth-order central difference scheme to improve the accuracy of spatial gradient calculation. The difference template is as follows:

[0160] ;

[0161] Deformation propagation velocity vector data In engineering applications, its maximum value along the roadway axis is usually further calculated. Location of peak occurrence , as a representative indicator of the transforming vanguard. When Exceeding the preset warning threshold When the speed is typically between 0.3 m / h and 0.5 m / h, the system marks the direction of the velocity vector as an arrow on the GIS visualization platform to help on-site managers intuitively judge the trend of deformation spread.

[0162] In another preferred embodiment of the present invention, the engineering parameters and weight calibration method for multi-source parameter normalized weighted mapping are as follows:

[0163] Design the limit geometric feature matrix The calibration was based on on-site roadway support design documents, measured reports of surrounding rock mechanical parameters, and historical disaster statistics, and was jointly calibrated by a group of technical experts in the mining area. Conservative values ​​(i.e., 10%–20% smaller) were used for each threshold value during the initial stage of system operation, and were subsequently adjusted using measured data over 3–6 months of system operation.

[0164] Inclination weight correction function The specific expression can be selected according to the characteristics of the fractures in the mining area, except as mentioned above. Besides the functional form, a piecewise linear form can also be used: when When in [0°, 30°] ;when When in [30°, 60°] Linear growth to 1.3; when When in [60°, 90°] The linear dropback to 1.0. This segmented form better reflects the engineering experience regarding the severity of shear cracks on-site.

[0165] Normalized weighted expert vector The initial calibration preferably adopts the Analytic Hierarchy Process (AHP). A group of technical experts from the mining area compares the importance of each of the four characteristic parameters pairwise, constructs a judgment matrix, and then solves for the eigenvector corresponding to the largest eigenvalue. This eigenvector is then normalized and used as the initial value of W. After the system is running, W is dynamically updated according to an adaptive online calibration mechanism.

[0166] Warning distance The preferred method is based on the magnitude of the deformation propagation velocity vector data. Dynamic adjustment: when hour, Take 5m; when hour, Take 10m; when When ≥1.0 m / h, Take 20m. Initial reference pumping pressure for flexible mold pump grouting filling equipment. Typically, the pressure is taken as 60% to 70% of the equipment's rated pressure, after reverse stiffness compensation. The pressure should not exceed 95% of the equipment's rated pressure to protect equipment safety. The pressure gain coefficient β should initially be a conservative value of β=0.3, and gradually adjusted to the optimal value based on the measured response characteristics of the filling material on-site.

[0167] Industrial fieldbus selection should prioritize protocols that support time synchronization and have strong electromagnetic interference resistance (such as PROFINET IRT or EtherCAT), and reserve at least 3 levels of redundant transmission paths to ensure reliable issuance of control commands in complex electromagnetic environments underground.

[0168] The following is a complete implementation example of applying the method of this invention to a high-extraction fully mechanized longwall face:

[0169] The monitoring point is the transport roadway of the 8.8m high fully mechanized longwall mining face in a certain mine. The roadway has a cross-sectional size of 6.0m×4.5m. 32 displacement sensing nodes (spaced 1.0m apart), 8 principal stress sensors, 3 sets of roof deep base point separation instruments, and 12 filling body displacement sensors are arranged along the roadway axis. Two sets of hydraulic fracturing boreholes (10m deep) are drilled in the roof 20m in front of the working face. 4-oriented distributed fiber optic grating bundles are laid in the boreholes.

[0170] The system connects to each sensor node in real time via a PROFINET network, with a timestamp synchronization accuracy better than 100μs. The derivative of the principal stress increment time-series data is obtained... A total of 14 instantaneous loading events caused by coal mining machine cutting were identified on that day, of which 9 were... Value exceeds preset threshold This process involves identifying and removing abrupt changes in the original displacement data. Based on the Burgers rheological constitutive model, a time-domain hardening integral is performed on the baseline elastic displacement sequence after abrupt changes are removed, outputting pure rheological creep linear deformation data. The maximum value was found at a measuring point 8 meters in front of the workpiece. =42mm.

[0171] Obtain the measured value of the total linear displacement of macroscopic subsidence on the top surface. =68mm, take the geometric coupling calibration coefficient. (Conservative approximation), then the local delamination size =42mm, based on macroscopic total subsidence measurement value and =42mm, decoupling yields the background rigid body settlement data of the strata. =26mm; based on the original linear deformation of the axial section of the filling material. 58mm yielded local net compressive deformation data =32mm; Cooperative compression ratio data =32 / 42=0.76, which falls within the warning range. (in =1.0), the system triggers an on-site inspection notification.

[0172] The system performs temperature compensation (daily hydration heat temperature increment ΔT = 8℃) on the wavelength drift data of the 4-oriented fiber grating bundle and then inputs the strain space coordinate transformation matrix. Obtain the normal width data of the crack. =3.2mm, crack dip angle determined by E3D principal strain plane normal orientation. =52°. Seam width weight correction factor. =1+0.3· 52° =1.186, seam width after weighted adjustment =3.79mm.

[0173] Based on pure rheological creep linear deformation data from 32 measurement points Construct a 32×600 spatiotemporal sequence matrix (Δt=1min, duration 10h), and obtain it through fourth-order central difference. and Substitute The deformation propagation velocity vector is obtained. =0.42m / h, the direction is in front of the working face (i.e., in the negative direction along the roadway axis).

[0174] System by normalized weights Calculate the multidimensional dimensionless geometrically degenerate eigenvector Dimensionless comprehensive deformation characteristic index was obtained. .current =1.0, the system determines that the surrounding rock-fill body bearing system is in a safe state, only triggers the early warning, and does not issue stiffness compensation control commands.

[0175] Subsequently, after a blasting disturbance in an adjacent section... The value rose to 1.18, exceeding the safety critical control threshold. The system extracted... Directional symbols are used to define the boundary of the hazardous space as the section 8m to 18m in front of the work area. The reverse stiffness compensation control command was issued to the flexible mold pumping grouting and filling equipment in this section, so that... From the initial reference pressure The pressure has been increased to 2.18 MPa.

[0176] After the reverse stiffness compensation control command is issued, the system continuously monitors the newly acquired displacement and stress increment data and calculates the expected rheological displacement. At the 3rd hour, the evolution rate of the actual deformation reconstruction residual r(t) was detected to be ṙ = 0.62 mm / h, exceeding the preset convergence threshold. =0.5mm / h, the system determines that the surrounding rock has generated a nonlinear accelerated failure disturbance, retrieves the current time Δσ(t)=0.4MPa as the disturbance correction factor κ=0.4 / 0.3=1.33, and applies it to the normalized weighted expert vector. +0.02×1.33×( The linear adaptive update of the shape variable weights The propagation speed weight has been increased from 0.35 to 0.38. The value was increased from 0.20 to 0.22, with other weights adjusted accordingly and normalized to make subsequent indicator calculations more consistent with the current accelerated instability characteristics.

[0177] This implementation example verifies the comprehensive effectiveness of the method of the present invention in monitoring the deformation of the surrounding rock in high-intensity mining roadways: the accuracy rate of identifying abrupt change components is not less than 95%; the decoupling error of the formation background rigid body sedimentation is controlled within ±3mm; the normal width measurement accuracy of the distributed fiber grating crack geometric inversion is better than ±0.2mm, and the tilt angle measurement accuracy is better than ±3°; the dimensionless comprehensive deformation characteristic index provides an early warning time of not less than 2 hours for actual disaster events; the system has been operating stably for more than 6 months without any false alarms or missed alarms caused by background rigid body displacement interference, providing reliable data support for mine pressure management and scheduling.

[0178] Example 2:

[0179] like Figure 3 and Figure 4 As shown, the coal mine roadway surrounding rock deformation monitoring system disclosed in this invention uses the aforementioned coal mine roadway surrounding rock deformation monitoring method. The system is integrated into the edge computing node of the mine's integrated automation platform or the ground monitoring center server. Through a three-tier architecture of distributed sensing terminals, underground edge computing units, and a ground analysis platform, it constructs a real-time surrounding rock stability monitoring system with multi-source data fusion and closed-loop control capabilities, specifically including the following modules:

[0180] The data stream acquisition module is used to acquire discretely distributed raw linear displacement time-series data, principal stress increment time-series data, reflected Bragg wavelength drift time-series data, macroscopic total subsidence linear displacement of the roof surface, and axial raw linear deformation of the filling body along the tunnel axis. This module establishes a connection with the data acquisition front-end composed of underground sensor terminals via an industrial Ethernet / PROFINET network. It achieves timestamp alignment of each sensor node based on the IEEE 1588 precision clock synchronization protocol and maintains a raw data buffer circular queue to support time-series backtracking calculations for subsequent modules. In terms of hardware implementation, this module relies on a multi-channel isolated data acquisition interface card of the edge computing node (supporting RS-485, 4-20mA current loop, and fiber optic communication) to ensure the reliability and intrinsic safety of data acquisition in the strong electromagnetic interference environment underground.

[0181] The stripping and solving module is used to remove abrupt components from the original linear displacement time series data based on the instantaneous loading rate of the principal stress increment time series data. Then, using this principal stress increment time series data as the integration stress path, it performs rheological constitutive integral solving to obtain pure rheological creep linear deformation data. This module incorporates a preset rate-of-change threshold online sliding update mechanism and a rheological constitutive compliance kernel function calibration parameter library. The time-domain hardening integral operation is executed in floating-point parallel mode on the multi-core ARM processor of the edge computing unit, with a single frame data processing time of less than 20ms, meeting the real-time requirements of surrounding rock deformation monitoring.

[0182] The calculation module uses pure rheological creep linear deformation data as a local shallow delamination reference. It performs spatial geometric difference calculations with the macroscopic total subsidence linear displacement to obtain the formation background rigid body settlement data. Then, it performs difference calculations with the axial original linear deformation data and the formation background rigid body settlement data to measure the local net compressive deformation data of the infill body, and finally calculates the cooperative compression ratio data. This module achieves physical decoupling between formation background rigid body settlement and local net compressive deformation, and maintains preset support design cooperative benchmark values. The rolling correction logic; it also integrates a three-level collaborative compression ratio determination submodule, based on collaborative compression ratio data. The numerical values ​​automatically trigger three status indicators: normal, warning, and danger, and upload the judgment results to the mining pressure management and scheduling system via the OPC UA protocol.

[0183] The crack geometry inversion module is used to input the reflected Bragg wavelength drift time series data into a preset strain space coordinate transformation matrix for inverse transformation calculation, obtaining crack normal width data and crack dip angle data. This module is based on the high-precision wavelength measurement results of a distributed fiber optic grating demodulator (wavelength scanning type), and calculates the crack width and dip angle based on the photoelectric strain sensitivity coefficient. The strain is mapped to a multi-orientation line strain vector, and then the three-dimensional spatial strain tensor is reconstructed through Moore-Penrose generalized inverse matrix operations. The module incorporates temperature compensation sub-units and multi-plane piecewise inversion sub-units, maintaining inversion accuracy under conditions of hydration heat temperature drift and complex crack morphology. (Strain space coordinate transformation matrix) The construction parameters and multi-plane segmentation trigger conditions are stored in the module's non-volatile memory, supporting remote configuration updates.

[0184] The solution module constructs a spatiotemporal sequence of surrounding rock geometric deformation data from pure rheological creep linear deformation data and performs differential gradient calculations to convert it into spatial geometric change rate and temporal rheological rate values. These values ​​are then substituted into the velocity calculation formula to obtain the deformation propagation velocity vector data. This module uses a fourth-order central difference operator as the core of the spatiotemporal grid gradient calculation, and solves for the deformation propagation velocity vector based on characteristic line theory. The calculation process is executed on a GPU or a dedicated SIMD instruction set acceleration unit, with an overall processing time of less than 100ms for a spatiotemporal sequence matrix of 32 measurement points × 600 time points. The solution results are uploaded to a ground GIS visualization platform via a web service interface, visually presenting the deformation front direction and propagation rate as arrows.

[0185] The index evaluation module extracts and normalizes weighted data from co-compression ratio, crack normal width, crack dip angle, deformation propagation velocity vector data, and pure rheological creep linear deformation data to obtain dimensionless comprehensive deformation characteristic indices. This module includes a built-in design limit geometric characteristic matrix. The parameter management subunit of the normalized weighted expert vector W supports initial calibration based on the AHP (Analytic Hierarchy Process) and sliding correction based on operational data. The index calculation results are transmitted in real time to the underground audible and visual alarm terminal via the Modbus-TCP protocol and output to the mine safety management cloud platform via REST API.

[0186] Furthermore, this system preferably also includes a physical closed-loop control module and an environmental adaptive online calibration module:

[0187] The physical closed-loop control module performs numerical comparison logic operations between the dimensionless comprehensive deformation characteristic index and the preset safety critical control threshold. This identifies deformation exceeding limits and locks the dangerous space boundary. It then sends a reverse stiffness compensation control command to the inverter drive of the flexible mold pumping grouting filling equipment deployed within the dangerous space boundary, increasing the initial pumping filling pressure. This module establishes a connection with the inverter of the flexible mold pumping grouting filling equipment via PROFINET IRT industrial real-time Ethernet, ensuring an end-to-end transmission delay of less than 10ms for control commands, guaranteeing a rapid response from monitoring data triggering to on-site execution. The module also includes a built-in warning distance. The dynamic adjustment and pressure gain coefficient β protection limit logic prevents control command oscillations caused by drastic fluctuations in indicators.

[0188] The environmental adaptive online calibration module continuously acquires newly collected original linear displacement and principal stress increment data within the hazardous space boundary after the issuance of the reverse stiffness compensation control command. It then calls a pre-calibrated surrounding rock creep constitutive model to calculate the expected rheological displacement and identifies nonlinear accelerated failure disturbances through the actual deformation reconstruction residuals and their evolution rate values. Finally, it performs linear adaptive gain compensation updates on the normalized weighted expert vector. This module achieves real-time correction of the W vector based on online matrix operations and writes the updated results back to the parameter management subunit of the index evaluation module, forming a complete adaptive closed loop of "measurement-discrimination-control-calibration".

[0189] This invention achieves high-precision multi-physics field monitoring and closed-loop control of coal mine roadway surrounding rock deformation through multi-source sensor fusion at the hardware level and multi-module collaboration at the software level. It eliminates the contamination of real creep information by instantaneous elastic jumps and high-frequency vibration noise, reliably separates overall large deformation from local net deformation under high-extraction mining scenarios, simultaneously outputs dual geometric parameters of crack normal cracking and shear dip angle, unifies the instability discrimination criteria for multiple heterogeneous monitoring parameters, and realizes a complete closed loop from monitoring, discrimination, control to self-calibration. The technical effects of each module are mutually enhanced, improving the engineering practicality and intrinsic safety level of the coal mine roadway surrounding rock deformation monitoring system under extreme conditions such as high-extraction mining.

[0190] The technical scope of this invention is not limited to the content described above. Those skilled in the art can make various modifications and variations to the above embodiments without departing from the technical concept of this invention, and all such modifications and variations should fall within the protection scope of this invention.

Claims

1. A method for monitoring the deformation of surrounding rock in coal mine roadways, characterized in that: Includes the following steps: Acquire the original linear displacement time series data, principal stress increment time series data, reflected Bragg wavelength drift time series data, macroscopic total subsidence linear displacement of the roof surface, and axial original linear deformation of the filling body, which are discretely distributed along the roadway axis. Based on the instantaneous loading rate obtained by first-order time-domain differentiation of the principal stress increment time series data, the step abrupt component corresponding to the abnormal loading rate in the original linear displacement time series data is identified and removed to obtain the reference elastic displacement time series sequence. Then, using the principal stress increment time series data as the integral stress path, the rheological constitutive integral solution is performed on the reference elastic displacement time series to obtain pure rheological creep linear deformation data; The pure rheological creep linear deformation data is used as a local shallow delamination reference. The calibrated linear geometric coupling relationship is mapped to local delamination size data. Then, spatial geometric difference operation is performed with the macroscopic total subsidence linear displacement to obtain the formation background rigid body settlement data. Then, the difference operation is performed between the axial original linear deformation and the formation background rigid body settlement data to measure the local net compression deformation data of the infill body. Finally, the cooperative compression ratio data is calculated. The reflected Bragg wavelength drift time series data is input into a preset strain space coordinate transformation matrix for inverse transformation calculation to obtain crack normal width data and crack dip angle data determined by the principal strain plane orientation of the three-dimensional spatial geometric strain tensor. The spatial-temporal sequence data of surrounding rock geometric deformation is constructed by the pure rheological creep linear deformation data and differential gradient operation is performed to convert it into spatial geometric change rate value and temporal rheological rate value. Then, it is substituted into the characteristic linear velocity solution formula based on the advection model to obtain deformation propagation velocity vector data. The cooperative compression ratio data, crack normal width data, crack dip angle data, deformation propagation velocity vector data, and pure rheological creep linear deformation data are extracted, normalized and weighted, and dimensionless comprehensive deformation characteristic indexes are obtained by measurement.

2. The method for monitoring the deformation of surrounding rock in coal mine roadways according to claim 1, characterized in that: Obtaining pure rheological creep linear deformation data specifically includes the following steps: Extract the principal stress increment time series data and perform a first-order time-domain derivative operation on it to obtain the instantaneous stress loading rate; Identify the abnormal instantaneous moment when the instantaneous stress loading rate exceeds a preset rate of change threshold, and lock the step jump deformation size corresponding to the abnormal instantaneous moment in the original linear displacement time series data, define it as the abrupt change component and remove it, and output the reference elastic displacement time series sequence; Using the principal stress increment time series data as the stress-rheological integration path, time-domain hardening integration is performed on the reference elastic displacement time series to eliminate high-frequency symmetrical mechanical vibration noise with zero mean, thereby generating the pure rheological creep linear deformation data used to characterize the true rheological response of the surrounding rock.

3. The method for monitoring deformation of surrounding rock in coal mine roadways according to claim 1, characterized in that: The local net compression deformation data of the filling body is obtained by measurement, and then the cooperative compression ratio data is calculated. The specific steps include: The pure rheological creep linear deformation data is used as a local shallow delamination reference and mapped to the local delamination size data of the shallow unstable rock layer deviating from the deep stable anchor point according to the calibrated linear geometric coupling relationship. The linear geometric coupling relationship is that the local delamination size data is equal to the product of the geometric coupling calibration coefficient and the pure rheological creep linear deformation data. Based on the local delamination size data, the background rigid body settlement data of the strata is obtained by subtraction and decoupling. The difference operation is performed between the original axial linear deformation and the background rigid body settlement data of the stratum to obtain the local net compression deformation data that characterizes the actual change in the pressure of the filling body itself. The local net compression deformation data is used as the numerator, and the local delamination size data is used as the denominator. A division ratio operation is performed to output the cooperative compression ratio data, which is used to quantitatively characterize the reverse support constraint state between the top plate surface and the filling body.

4. The method for monitoring deformation of surrounding rock in coal mine roadways according to claim 1, characterized in that: Obtaining crack normal width and crack dip angle data involves the following steps: The reflected Bragg wavelength drift time series data is linearly mapped by a preset photoelectric strain sensitivity coefficient to obtain a one-dimensional multi-orientation line strain measurement sequence along the fiber axis. The one-dimensional multi-orientation line strain measurement sequence is input into the strain space coordinate transformation matrix to perform inverse transformation matrix multiplication operation, and the three-dimensional spatial geometric strain tensor measurement data is inverted and reconstructed; the strain space coordinate transformation matrix is ​​constructed by calling up the solid mechanics geometric compatibility equation based on the spatial orthogonal or cross-orientation geometric topology of the distributed fiber grating bundles laid in the hydraulic fracturing borehole. By introducing the discontinuous boundary condition of crack displacement in fracture mechanics, the opening-type geometric discontinuity displacement and the shear-type geometric discontinuity displacement are separated and measured from the three-dimensional spatial geometric strain tensor measurement data. These are respectively quantified and calculated as the crack normal width measurement data, which characterizes the linear size of crack opening, and the crack dip angle measurement data, which characterizes the crack surface relative to the axis of the measured surrounding rock, determined by the normal orientation of the principal strain plane of the three-dimensional spatial geometric strain tensor. Furthermore, the crack displacement mixing degree angle data, which characterizes the degree of transition from the opening-type to the shear-type geometric discontinuity displacement, is calculated from the arctangent of the ratio of the shear-type geometric discontinuity displacement to the opening-type geometric discontinuity displacement.

5. The method for monitoring the deformation of surrounding rock in coal mine roadways according to claim 4, characterized in that: Obtaining deformation propagation velocity vector data specifically includes the following steps: The pure rheological creep linear deformation data generated at multiple measuring points discretely distributed along the tunnel axis are spatially aligned, and a unified high-precision clock stamp is introduced to construct and generate the spatiotemporal sequence data of the surrounding rock geometric deformation containing two-dimensional geometric features of spatial axis and time axis. The preset spatiotemporal grid difference operator is invoked as the execution carrier for the difference gradient operation. The adjacent discrete grid difference calculation along the roadway spatial axis is performed on the spatiotemporal sequence data of the surrounding rock geometric deformation, and the spatial geometric change rate value is output. Simultaneously, the temporal adjacent time step difference calculation along the time axis is performed, and the temporal rheological rate value is output. The pure rheological creep deformation field of the surrounding rock is approximated as a kinematic wave transported unidirectionally along the tunnel axis. The characteristic line theory of mathematical physics equations is introduced as the velocity calculation formula. The derivative ratio of the time rheological rate value and the spatial geometric change rate value is negatively mapped to obtain the deformation propagation velocity vector data with spatial direction indication characteristics, which is used to characterize the geometric characteristics of rock mass deformation propagation. The axial vector direction of the deformation propagation velocity vector data is defined according to the positive or negative sign of the spatial geometric change rate value.

6. The method for monitoring deformation of surrounding rock in coal mine roadways according to claim 5, characterized in that: Before constructing the spatiotemporal sequence data of surrounding rock geometric deformation, a step to enhance the spatial consistency of deformation data is also included, specifically: The spatial coordinates of the total macroscopic subsidence linear displacement of the roof surface and the original axial linear deformation of the filling body are obtained. The spatial coordinates are determined and recorded by a total station with the axial mileage of the roadway as the unified spatial reference coordinates. Based on the sampling position of the pure rheological creep linear deformation data, spatial resolution interpolation resampling is performed on the total macroscopic subsidence linear displacement and the original axial linear deformation to achieve grid dimension alignment of multi-source heterogeneous deformation data. Simultaneously collect local temperature increment data caused by the heat release of underground filling hydration, and perform temperature compensation calculation on the one-dimensional multi-orientation line strain sequence based on the preset optical fiber thermal expansion correction coefficient, output the temperature-compensated one-dimensional multi-orientation line strain sequence, and replace the original one-dimensional multi-orientation line strain sequence in the strain space coordinate transformation matrix.

7. The method for monitoring the deformation of surrounding rock in coal mine roadways according to claim 1, characterized in that: The measurement of dimensionless comprehensive deformation characteristic index includes the following steps: Retrieve a preset design limit geometric feature matrix, which includes a limit deformation threshold, a limit compression ratio threshold, a limit propagation speed threshold, and a limit slit width threshold; Obtain a preset dip angle weight correction function, input the crack dip angle data into the dip angle weight correction function, and calculate the crack width weight correction coefficient; Multiply the crack normal width data by the crack width weight correction coefficient to obtain the weighted corrected crack normal width data; The modulus values ​​in the pure rheological creep linear deformation data, the cooperative compression ratio data, the deformation propagation velocity vector data, and the weighted corrected crack normal width data are respectively subjected to matrix division normalization operations with the corresponding limit deformation threshold, limit compression ratio threshold, limit propagation velocity threshold, and limit crack width threshold to generate a multidimensional dimensionless geometric degradation feature vector. Extract the preset normalized weighted expert vector, perform matrix inner product multiplication and addition mapping between the multidimensional dimensionless geometric degradation feature vector and the normalized weighted expert vector, and output the dimensionless comprehensive deformation feature index.

8. The method for monitoring deformation of surrounding rock in coal mine roadways according to claim 7, characterized in that: After obtaining the dimensionless comprehensive deformation characteristic index by measurement, the process also includes a physical closed-loop control command issuance and execution step, specifically including: Perform a numerical comparison logic operation between the dimensionless comprehensive deformation characteristic index and the preset safety critical control threshold to identify deformation exceeding events where the dimensionless comprehensive deformation characteristic index is greater than the safety critical control threshold. In response to the deformation exceeding the limit event, the spatial direction symbol contained in the deformation propagation velocity vector data is extracted, and the dangerous spatial boundary pointed to by the deformation propagation front in the roadway space is locked. Based on the difference between the dimensionless comprehensive deformation characteristic index and the safety critical control threshold, a reverse stiffness compensation control command is issued to trigger secondary reinforcement.

9. The method for monitoring the deformation of surrounding rock in coal mine roadways according to claim 8, characterized in that: It also includes an environment-adaptive online calibration step, specifically including: After issuing the reverse stiffness compensation control command, the original linear displacement time series data and principal stress increment time series data re-acquired within the dangerous space boundary are continuously acquired. Based on the principal stress increment time series data at the current moment, the pre-calibrated surrounding rock creep constitutive model and its parameter set are called to calculate the actual stress increment change rate after reverse stiffness compensation, and the data is substituted into the surrounding rock creep constitutive model for time domain integration to generate the expected rheological displacement. Calculate the measured displacement increment of the original linear displacement time series data at the current time step, and calculate the actual deformation reconstruction residual between it and the expected rheological displacement. The actual deformation reconstruction residual is subjected to a first-order difference operation with respect to the time dimension to obtain the residual evolution rate value; when the residual evolution rate value is greater than the preset convergence threshold, it is determined that the surrounding rock has generated a nonlinear accelerated failure disturbance, the principal stress increment time series data at the current moment is retrieved as the disturbance correction factor, and the normalized weight expert vector is updated with linear adaptive gain compensation.

10. A coal mine roadway surrounding rock deformation monitoring system, characterized in that: The method for monitoring the deformation of surrounding rock in coal mine roadways as described in any one of claims 1 to 9 includes: The data stream acquisition module is used to acquire the original linear displacement time series data, principal stress increment time series data, reflected Bragg wavelength drift time series data, macroscopic total subsidence linear displacement of the roof surface, and axial original linear deformation of the filling body, which are discretely distributed along the tunnel axis. The stripping and solving module is used to identify and remove step abrupt components corresponding to abnormal loading rate moments in the original linear displacement time series data based on the instantaneous loading rate obtained by first-order time-domain differentiation of the principal stress increment time series data, thereby obtaining a reference elastic displacement time series sequence; then, using the principal stress increment time series data as the integral stress path, the reference elastic displacement time series sequence is subjected to rheological constitutive integral solving to obtain pure rheological creep linear deformation data. The calculation module is used to use the pure rheological creep linear deformation data as a local shallow delamination reference, map it to local delamination size data through the calibrated linear geometric coupling relationship, and then perform spatial geometric difference operation with the macroscopic total subsidence linear displacement to obtain the formation background rigid body settlement data. Then, the original axial linear deformation variable is performed with the formation background rigid body settlement data to measure the local net compression deformation data of the infill body, and then the cooperative compression ratio data is calculated. The crack geometry inversion module is used to input the reflected Bragg wavelength drift time series data into a preset strain space coordinate transformation matrix for inverse transformation calculation to obtain crack normal width data and crack dip angle data determined by the principal strain plane orientation of the three-dimensional spatial geometric strain tensor. The solution module is used to construct the spatiotemporal sequence data of surrounding rock geometric deformation through the pure rheological creep linear deformation data and perform differential gradient calculation to convert it into spatial geometric change rate value and temporal rheological rate value, and substitute it into the characteristic linear velocity solution formula based on the advection model to obtain deformation propagation velocity vector data. The index evaluation module is used to extract the cooperative compression ratio data, crack normal width data, crack dip angle data, deformation propagation velocity vector data, and pure rheological creep linear deformation data, perform normalized weighted mapping, and measure to obtain a dimensionless comprehensive deformation characteristic index.

Images