Coal mine area goaf geological disaster forecasting and early warning system and method

Abstract

The invention belongs to the technical field of coal mining, and discloses a coal mine area goaf geological disaster forecasting and early warning system and method. The method comprises the following steps: deploying a ground surface monitoring layer to calculate a ground surface settlement gradient field in real time; deploying an underground sensing layer to collect rock mass stress-acoustic emission signals in real time, and collecting gas concentration data; a dynamic fracture network evolution model is used for obtaining the fractal dimension of the fracture in the goaf, the fracture extension trend is inversed through microseismic data, and the seepage-stress relation of the goaf is obtained; and inputting the obtained information into a bee colony optimization neural network (BEO-NN), and obtaining graded early warning information of geological disasters in the goaf of the coal mine area through a bionic intelligent early warning mechanism. According to the method, the coal mine geological disaster hidden danger can be more effectively identified and pre-warned, so that the safety and sustainability of coal mine production are guaranteed.

Description

Technical Field

[0001] This invention belongs to the field of coal mining technology, and in particular relates to a geological disaster forecasting and early warning system and method for goaf areas in coal mines. Background Technology

[0002] Coal mining is primarily underground, and the process disrupts the balance of the original geological formations, leading to geological disasters. These disasters not only cause economic losses to coal mining companies but also result in significant casualties and have extremely negative social impacts. The cause of these disasters is that mining disrupts the balance of the original geological formations, and abnormal phenomena occur during the process of achieving a new balance. Generally, coal mine geological disasters can be categorized into those caused by mining techniques and those caused by uncertainties in geological conditions. To prevent coal mine geological disasters, efforts should be made to conduct thorough coal mine geological exploration, adopt reasonable mining methods, and focus on gas extraction from coal seams.

[0003] Existing patented technologies cannot simultaneously achieve comprehensive detection, early warning, and information management of typical natural disasters such as mine disasters, geological disasters, and environmental disasters. This makes it impossible to realize dynamic monitoring and real-time forecasting and early warning of disaster risks and hazards, resulting in poor disaster detection, early risk identification, and forecasting and early warning capabilities. Improving the uniformity, accuracy, and reliability of disaster data acquisition and management is an important direction for current research and development.

[0004] In the coal mining sector, geological hazards are a significant factor restricting safe production and sustainable development in coal mines. Due to the complexity of the geological environment in coal mines and the constant disturbance caused by mining activities, geological hazards such as ground subsidence seriously threaten the lives and property of coal mine workers and the normal production and operation of coal mines.

[0005] Traditionally, the monitoring and early warning of geological disasters in coal mines have relied primarily on manual inspections and simple monitoring equipment, a method with numerous limitations. First, manual inspections struggle to achieve comprehensive, real-time monitoring of the coal mine geological environment, leading to frequent missed and false detections. Second, simple monitoring equipment often provides only single geological indicator data, failing to comprehensively assess the risk and development trend of geological disasters. Therefore, traditional methods are significantly inadequate in terms of the accuracy and timeliness of geological disaster early warning.

[0006] In recent years, with the rapid development of Internet of Things (IoT) technology, its applications in various fields have become increasingly widespread. IoT technology, through intelligent sensing, identification, and ubiquitous computing, forms a vast network that combines various information sensing devices with the internet. This technology provides a new solution for the monitoring and early warning of geological disasters in coal mines.

[0007] Furthermore, traditional monitoring relies on single methods (such as InSAR surface deformation monitoring and gas concentration analysis), which fails to achieve data collaboration between surface and downhole operations, resulting in delayed early warnings; sensors are susceptible to interference from high humidity and high dust levels, leading to high data distortion rates; and early warning models are static and do not consider the dynamic coupling effect between mining stress and geological structure. Summary of the Invention

[0008] To overcome the problems existing in related technologies, the present invention discloses an embodiment of a geological disaster forecasting and early warning system and method for goaf areas in coal mines.

[0009] The technical solution is as follows: A method for forecasting and early warning of geological disasters in coal mining goaf areas, comprising the following steps:

[0010] S1, deploys an aboveground monitoring layer to calculate the surface subsidence gradient field in real time, and deploys a downhole sensing layer to collect rock stress-acoustic emission signals in real time. Gas concentration data;

[0011] S2. The fractal dimension of the fractures in the goaf is obtained by using a dynamic fracture network evolution model, and the fracture propagation trend is inverted by microseismic data to obtain the relationship between seepage and stress in the goaf.

[0012] S3 inputs the acquired information into the bee colony optimization neural network BSO-NN, and obtains graded early warning information on geological disasters in the goaf of the coal mining area through a biomimetic intelligent early warning mechanism.

[0013] In step S1, the improved millimeter-wave interferometric radar operates at a frequency of 80 GHz in the well monitoring layer.

[0014] The real-time calculation of the surface subsidence gradient field by integrating BeiDou / GNSS data includes: introducing deep learning YOLOv10 to automatically identify subsidence funnel regions, combining a two-dimensional Gaussian model to invert nonlinear deformation, and solving phase aliasing in large gradient subsidence; as well as the calculation of the surface subsidence gradient field.

[0015] Furthermore, the construction of the two-dimensional Gaussian model includes:

[0016] For the inversion of settling funnel parameters, a two-dimensional Gaussian model is constructed for each funnel region identified by YOLOv10, with the objective function of minimizing the angle deviation. The expression is as follows:

[0017]

[0018] In the formula, The area of ​​the funnel region. It is horizontal. Vertically, As the central location, To control the range of deformation;

[0019] The deformation gradient is inverted by fitting the parameters of a two-dimensional Gaussian model through nonlinear optimization.

[0020] Furthermore, the calculation of the surface subsidence gradient field includes:

[0021] Spatiotemporal interpolation and gridding: The fused displacement data is discretized according to a spatiotemporal grid, and Kriging interpolation or inverse distance weighting (IDW) is used to fill the gaps.

[0022] Gradient field generation includes: calculating settlement rate based on time-series displacement data, per grid cell, with the expression:

[0023]

[0024] In the formula, For the time difference, It is the displacement difference;

[0025] Gradient vector extraction involves calculating the spatial gradient using the Sobel operator or the central difference method, expressed as:

[0026]

[0027] Characterizes the intensity and direction of spatial variation in settling rate.

[0028] In step S1, the downhole sensing layer includes: a biomimetic anchor bolt group sensor, which simulates the structure of plant roots and is embedded with a piezoelectric ceramic and fiber optic grating composite sensing unit to collect rock stress-acoustic emission signals in real time;

[0029] The biomimetic anchor bolt group sensor includes:

[0030] Strain-acoustic emission dual-mode sensing: fiber optic gratings or vibrating wire micro strain gauges are embedded inside the anchor bolt to monitor stress changes, and piezoelectric ceramic acoustic emission sensors are integrated on the surface of the bolt to capture high-frequency elastic waves from rock fracture.

[0031] The grid-like layout, with multiple sets of sensors distributed along the anchor bolt axis, forms a spatial sensing network to achieve three-dimensional mapping of the stress field.

[0032] The biomimetic structural design references the resonance principle of biological cilia to design a biomimetic anchor group sensor.

[0033] Real-time acquisition of rock mass stress-acoustic emission signals includes:

[0034] (1) Signal synchronization and noise reduction: Time synchronization of group sensors is achieved through GPS or high-precision crystal oscillator; wavelet transform and EMD decomposition are used to eliminate interference from downhole electromechanical equipment while retaining effective acoustic emission characteristics;

[0035] (2) Edge intelligent processing: deploy a lightweight LSTM network on the sensor nodes of the bionic anchor group to identify acoustic emission event types and associate stress mutation points in real time; dynamically compress data and upload only the characteristic parameters of abnormal events;

[0036] (3) Multimodal data fusion; establish a stress-acoustic emission correlation model, and issue a warning of local overload of rock mass when stress suddenly increases and high-frequency acoustic emission events are added; when stress gradually changes and continuous low-energy acoustic emission is added, creep damage accumulation is detected.

[0037] The real-time acquisition of rock mass stress-acoustic emission signals is also corrected in real time using the Mie scattering model, including:

[0038] Input parameters: dust concentration, particle size distribution; calculation process:

[0039]

[0040] In the formula, Historical dust scattering concentration This represents the actual measured dust scattering concentration. The dust scattering coefficient is... For optical path, This refers to the dust concentration.

[0041] In step S2, the dynamic fracture network evolution model includes:

[0042] Establish a model for calculating the fractal dimension of fractures in goaf areas:

[0043]

[0044] In the formula, For crack scanning accuracy, To determine the number of fractures, the fracture propagation trend is inverted using microseismic data.

[0045] The coupled seepage-stress equation is:

[0046]

[0047] In the formula, For penetration rate, penetration rate With stress Related, Pore ​​pressure, The saturation level is used to predict the risk of water inrush / gas outburst. For the vector differential operator of the downward gradient, Porosity For a specific moment;

[0048] Inversion of fracture propagation trends using microseismic data includes:

[0049] S201, Data Acquisition and System Deployment: A three-component geophone array is deployed in the target area to form a microseismic monitoring network; the system acquires P-wave and S-wave elastic wave signals generated by rock fracture in real time, recording the time, amplitude and polarity.

[0050] S202, Data preprocessing; bandpass filtering of the raw seismic record; static and dynamic corrections;

[0051] S203, Microseismic event location; using the P-wave first arrival times of at least 4 sensors, the source coordinates are calculated using travel-time inversion or grid search methods;

[0052] S204, focal mechanism and parameter inversion; waveform type is determined by polarization analysis, and focal mechanism parameters, including strike, dip, and slip angle, are extracted; a multi-scan combined waveform inversion algorithm is used to simultaneously optimize spatial location and focal mechanism, and identify shear rupture direction; microseismic energy, apparent stress, and volumetric potential index are calculated to reflect rupture intensity and extent.

[0053] Furthermore, polarization analysis determines waveform type, including:

[0054] Three-component data matrix construction: A data matrix is ​​constructed from the three component data segments recorded by the seismograph: vertical Z, north-south N, and east-west E.

[0055]

[0056] Covariance matrix calculation:

[0057]

[0058] In the formula, superscript Let be the transpose matrix, and be the number of sampling points within the time window;

[0059] Eigenvalue decomposition: Solve for the eigenvalues ​​λ1≥λ2≥λ3 and the corresponding eigenvectors v1,v2,v3 of the covariance matrix; where, for P-wave determination, the eigenvalues ​​satisfy λ1≫λ2≈λ3;

[0060] The polarization direction, the eigenvector v1 corresponding to the maximum eigenvalue λ1 indicates the wave propagation direction; among them, for S-wave determination, the eigenvalues ​​satisfy λ1≈λ2>λ3;

[0061] The polarization complexity is high, and combined with the axial ratio of the polarization ellipse, the axial ratio = λ2 / λ1.

[0062] Furthermore, the focal mechanism parameters were obtained through inversion of the P-wave initial polarity or the entire waveform at multiple stations;

[0063] The P-wave initial polarity method includes:

[0064] Identify the P-wave initial motion direction; mark the initial motion symbols of each station in the lower hemisphere projection of the source sphere; calculate parameters by finding two orthogonal nodal planes so that the initial motion symbols are located in the four quadrants; the parameters are defined geometrically by the nodal planes.

[0065] The strike φ, the azimuth of the line of intersection of the fault and the horizontal plane relative to true north, ranging from [0°, 360°);

[0066] Dip angle δ, the angle between the fault plane and the horizontal plane, ranging from 0° to 90°;

[0067] Sliding angle λ is the angle between the sliding direction and the direction of travel, ranging from [-180°, 180°] to [0°, 360°].

[0068]

[0069] Wherein, λ≈0°, left-lateral strike-slip; λ≈±180°, right-lateral strike-slip; λ≈+90°, reverse fault; λ≈−90°, normal fault;

[0070] The full waveform inversion method includes:

[0071] Minimize the objective function using observed waveforms Compared with theoretical waveforms The residuals are solved by grid search or iterative optimization:

[0072]

[0073] In the formula, For model parameters, ;

[0074] P-wave and surface wave bands were separated and inverted separately using the CAP method, followed by weighted fitting:

[0075]

[0076] In the formula, The weighted fitting values ​​are obtained after inversion of the P-wave and surface wave bands respectively. All are weighting coefficients. To separate the observed waveform residuals from P-wave inversion, To separate the theoretical waveform residuals from P-wave inversion, For the observed waveform residuals of the separated surface wave inversion, The theoretical waveform residual for the separation of surface waves inversion;

[0077] Micro-vibration energy The calculation formula is based on the vibration velocity decay model:

[0078]

[0079] In the formula, The peak vibration velocity at the measuring point is obtained directly through the sensor. The distance from the seismic source to the sensor. This is an energy characteristic coefficient, which is related to the characteristics of the site medium.

[0080] Seismic Moment The physical quantity characterizing the fault slip intensity is calculated based on waveform inversion:

[0081]

[0082] In the formula, For the density of the medium, For seismic wave velocity, The low-frequency level of the seismic source spectrum was obtained through spectral analysis. For radiation mode coefficients, This is a directional correction factor;

[0083] Related indicators: seismic moment and moment magnitude The relationship is:

[0084]

[0085] Apparent stress The average stress consumed per unit area of ​​fault slip:

[0086]

[0087] Body change The expression representing the volumetric change caused by fault slip is:

[0088]

[0089] In the formula, This is the shear modulus.

[0090] In step S3, the bee colony optimization neural network BSO-NN includes:

[0091] Input layer: multidimensional heterogeneous data on surface deformation, gas concentration, microseismic energy, and fracture density;

[0092] Hidden layer, dynamic weight adjustment module, simulates the process of bee foraging path optimization;

[0093] The implementation process of risk classification and early warning includes:

[0094] Phase 1: Data acquisition and feature extraction; fusion of multi-source monitoring data, integrating multi-dimensional heterogeneous data such as surface deformation, gas concentration, microseismic energy, and fracture density to generate time series features; screening of key risk factors, using Pearson correlation coefficient analysis to analyze the correlation between indicators and disasters, and eliminating redundant parameters;

[0095] Phase 2: Construction of dynamic risk classification model; initial weight allocation; initial weights are generated using Tent chaotic mapping based on historical disaster data, and classification thresholds are set.

[0096] Phase 3: Real-time early warning and response; early warning triggering mechanism, which triggers multi-level response when the model outputs a risk level that jumps; dynamic threshold adjustment, which recursively updates weights based on new monitoring data to achieve model self-calibration.

[0097] Another objective of this invention is to provide a geological disaster forecasting and early warning system for goaf areas in coal mines. This system implements the aforementioned geological disaster forecasting and early warning method for goaf areas in coal mines. The system includes:

[0098] The data acquisition module is used to deploy an aboveground monitoring layer to calculate the surface subsidence gradient field in real time; and to deploy a downhole sensing layer to acquire rock mass stress-acoustic emission signals in real time, as well as to acquire... Gas concentration data;

[0099] The goaf seepage-stress relationship acquisition module is used to obtain the fractal dimension of the goaf fracture using a dynamic fracture network evolution model, invert the fracture propagation trend through microseismic data, and obtain the goaf seepage-stress relationship.

[0100] The graded early warning module is used to input the acquired information into the bee colony optimization neural network BSO-NN, and obtain graded early warning information on geological disasters in the goaf of coal mining areas through a biomimetic intelligent early warning mechanism.

[0101] Combining all the above technical solutions, the beneficial effects of this invention are as follows: by deploying various geological disaster data acquisition devices and using Internet of Things technology to achieve real-time monitoring and transmission of data, combined with advanced data processing and analysis technologies, this invention can more effectively identify and warn of potential geological disasters in coal mines, thereby ensuring the safety and sustainability of coal mine production. Attached Figure Description

[0102] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this disclosure and, together with the description, serve to explain the principles of this disclosure;

[0103] Figure 1 This is a flowchart of a method for forecasting and early warning of geological disasters in coal mining subsidence areas provided in an embodiment of the present invention;

[0104] Figure 2 This is a flowchart of the process of inverting fracture propagation trends using microseismic data, provided in an embodiment of the present invention. Detailed Implementation

[0105] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0106] Example 1, such as Figure 1 Methods for forecasting and early warning of geological disasters in coal mining subsidence areas include:

[0107] S1, deploys an aboveground monitoring layer to calculate the surface subsidence gradient field in real time, and deploys a downhole sensing layer to collect rock stress-acoustic emission signals in real time. Gas concentration data;

[0108] S2. The fractal dimension of the fractures in the goaf is obtained by using a dynamic fracture network evolution model, and the fracture propagation trend is inverted by microseismic data to obtain the relationship between seepage and stress in the goaf.

[0109] S3 inputs the acquired information into the bee colony optimization neural network BSO-NN, and obtains graded early warning information on geological disasters in the goaf of the coal mining area through a biomimetic intelligent early warning mechanism.

[0110] This invention provides a geological disaster forecasting and early warning system for goaf areas in coal mines, the system comprising:

[0111] The data acquisition module is used to deploy an aboveground monitoring layer to calculate the surface subsidence gradient field in real time; and to deploy a downhole sensing layer to acquire rock mass stress-acoustic emission signals in real time, as well as to acquire... Gas concentration data;

[0112] The goaf seepage-stress relationship acquisition module is used to obtain the fractal dimension of the goaf fracture using a dynamic fracture network evolution model, invert the fracture propagation trend through microseismic data, and obtain the goaf seepage-stress relationship.

[0113] The graded early warning module is used to input the acquired information into the bee colony optimization neural network BSO-NN, and obtain graded early warning information on geological disasters in the goaf of coal mining areas through a biomimetic intelligent early warning mechanism.

[0114] For example, the present invention integrates surface millimeter-wave interferometric radar, downhole biomimetic anchor bolt sensing network, dynamic fracture network evolution model and bee colony optimized neural network to achieve dynamic coupling early warning.

[0115] For example, the surface of the biomimetic anchor sensor is covered with a pyrophyllite-like dust-repellent coating, and the operating temperature range is -40℃ to 150℃.

[0116] The input data for the fractal dimension calculation model of the goaf fracture comes from the three-dimensional location results of microseismic events, with a location error ≤3m;

[0117] For example, the present invention constructs a four-dimensional perception network consisting of ground, air, and well; wherein, in step S1, the well monitoring layer: deploys an improved millimeter-wave interferometric radar (operating frequency band 80GHz), integrates BeiDou / GNSS data, and calculates the surface subsidence gradient field in real time with an accuracy of ±0.5mm;

[0118] For example, 80GHz high-frequency millimeter-wave radar technology; ultra-high frequency band and wide bandwidth (80GHz, 7GHz bandwidth); provides sub-millimeter ranging accuracy (wavelength only 3.75mm), significantly improving deformation monitoring resolution; 7GHz ultra-wide bandwidth enhances anti-interference capability, reduces multipath effects, and adapts to complex terrain environments.

[0119] An adaptive temperature compensation mechanism is used to correct hardware drift in real time through a built-in temperature sensor, reducing the Pearson correlation between the intermediate frequency signal amplitude and temperature drift by 84%; ensuring measurement stability in extreme environments and avoiding the accuracy degradation caused by temperature changes in traditional radar.

[0120] An example is a deep fusion architecture for multi-source data;

[0121] i) Spatiotemporal synchronization optimization (BeiDou time synchronization + radar-camera calibration): BeiDou nanosecond-level time synchronization is used to achieve time synchronization between radar and GNSS data, with a timing alignment error of <1ms; combined with camera calibration parameters, the spatial coordinate transformation error is controlled within ±5cm, improving the consistency of three-dimensional space.

[0122] ii) Sensor-level pre-fusion algorithm: The radar point cloud and GNSS location information are fused at the raw data layer, and global correlation features are extracted through the Transformer backbone network to solve the limitation of single sensor perspective; Compared with traditional post-fusion, the target association error rate is reduced by 23.4%.

[0123] For example, innovations in gradient field solution algorithms;

[0124] (a) Dynamic phase unwrapping and Gaussian inversion model; Deep learning (YOLOv10) is introduced to automatically identify the settlement funnel region, and a two-dimensional Gaussian model is combined to invert nonlinear deformation, solving the phase aliasing problem of large gradient settlement, and reducing the root mean square error by 57.4%;

[0125] For example, the specific steps include:

[0126] Step 1: Data preparation and sedimentation funnel identification;

[0127] Step 1.1, Multi-source dataset construction: Collect interferograms of simulated and real mining areas (such as Sentinel-1A and LuTan-1 data), covering typical subsidence funnel regions (such as Datong, Shanxi, and mining areas in Inner Mongolia). Annotate subsidence funnel regions: Use rectangular boxes to annotate the funnel locations in the interferograms, forming a YOLO format training set (category x_center y_center width_height). Data augmentation: Employ strategies such as brightness adjustment, Gaussian noise injection, and random flipping to improve model generalization.

[0128] Step 1.2, YOLOv10 model training and optimization; Network improvements: The DWR-C2f module is introduced to enhance multi-scale feature extraction capabilities and improve the detection accuracy of small targets (small sedimentation funnels). Anchor-Free architecture is adopted to reduce parameter redundancy, and dynamic label allocation optimizes positive and negative sample matching. Loss function: Inner-WIoU regression loss is used to alleviate the problem of sample difficulty imbalance and improve localization accuracy (experiments show mAP50 reaches 92.00%). Lightweight deployment: FP16 quantization is performed through TensorRT to compress the model size (102MB → 53MB) and improve inference speed;

[0129] Step 2, Two-dimensional Gaussian inversion modeling;

[0130] 2.1) Settlement funnel parameter inversion;

[0131] For each funnel region identified by YOLOv10, a two-dimensional Gaussian model is constructed with the objective function of minimizing the angular deviation:

[0132]

[0133] In the formula, The area of ​​the funnel region. It is horizontal. Vertically, As the central location, To control the range of deformation;

[0134] The deformation gradient is inverted by fitting the parameters of a two-dimensional Gaussian model through nonlinear optimization (such as the Levenberg-Marquardt algorithm).

[0135] 2.2) Model Phase Generation and Separation: The Gaussian model is converted into a phase field to generate the simulated settling funnel phase components. The model phase is subtracted from the original interferogram to obtain the residual phase field, which significantly reduces the local phase gradient.

[0136] Step 3, Phase optimization and untangling;

[0137] 3.1) Residual Phase Filtering: Adaptive filtering (such as Goldstein filtering) is applied to the residual phase to reduce phase aliasing caused by noise and improve the signal-to-noise ratio of the interferogram. The filtering strength needs to be matched with the local coherence to avoid excessive smoothing that leads to loss of detail.

[0138] 3.2) Phase unwrapping and recovery: The filtered residual phase is unwrapped using a traditional unwrapping algorithm (such as minimum cost flow). The unwrapped residual phase is then superimposed with the model phase to recover the complete true deformation phase.

[0139] Step 4, Validation and Performance Evaluation;

[0140] 4.11) Comparison of quantitative indicators; verification through simulation data: statistical analysis of the number of phase residuals and root mean square error (RMSE) before and after unwrapping. Practical applications (such as the Datong mining area) show that the average RMSE is reduced by 57.4% compared to the traditional method; the number of residuals is reduced by 23.4%.

[0141] 4.2) Optimization of actual deployment; Multi-source data fusion: Combining LiDAR and hyperspectral data to improve generalization ability in complex terrain. Edge computing: Deploying models on embedded platforms such as Jetson to achieve real-time monitoring (inference latency <23ms). A summary of technical advantages is shown in Table 1.

[0142] Table 1 Technical Advantages

[0143]

[0144] This method has been applied to subsidence monitoring in mining areas, successfully identifying 3,842 subsidence funnels and providing high-precision data support for geological disaster prevention and control. In the future, it can be extended to monitoring multiple scenarios such as landslides and urban subsidence. For example, this invention supports complementary verification using millimeter-wave radar point clouds and InSAR data, improving the reliability of monitoring complex terrain.

[0145] (b) Real-time Kalman filter fusion engine; fuses radar instantaneous displacement data with GNSS long-term trajectory, optimizes convergence speed through adaptive step size strategy, and improves solution efficiency by 40%; outputs surface subsidence gradient field vector map (including horizontal displacement and vertical velocity), supporting geological hazard risk assessment.

[0146] For example, fusing radar instantaneous displacement data with GNSS long-term trajectories and optimizing the convergence speed through an adaptive step-size strategy includes:

[0147] (b.1) Multi-source data synchronization and preprocessing;

[0148] b.1-1) Unified Spatiotemporal Reference; Time Synchronization: Using hardware PPS signals or software interpolation methods (such as linear / Lagrange interpolation), the timestamps of radar (sampling rate ≥100Hz) and GNSS (1-50Hz) data are aligned to millisecond-level accuracy to eliminate time delay errors. Coordinate Transformation: The radar local coordinate system data is transformed to the CGCS2000 or WGS84 geodetic coordinate system used by GNSS through rotation and translation matrices to ensure spatial consistency.

[0149] b.1-2) Data quality control; outlier filtering: remove transient noise (such as equipment vibration interference) from radar data based on the 3σ criterion or Hampel filter. GNSS drift correction: correct atmospheric delay error using base station differential data (RTK / PPP), with static monitoring accuracy required to reach ±1mm.

[0150] (b.2) Design and implementation of adaptive fusion algorithm;

[0151] (b.2-1) Dual-channel joint optimization framework;

[0152] b.2-1.1) Delay-displacement decoupling modeling includes: Delay tracking unit: Based on the long-term GNSS trajectory, the system delay deviation D between the radar and GNSS is dynamically estimated using the LMS algorithm.

[0153]

[0154]

[0155] In the formula, For the first Node system delay deviation The residuals under systematic bias, For the delay of radar equipment;

[0156] Displacement fusion unit: The RLS algorithm is used to adaptively compensate for sensor gain differences and output fused displacement Sfused.

[0157] b.2-1.2) Variable step size strategy accelerates convergence; step size is dynamically adjusted based on the residual rate of change. Nonlinear adjustment step size :

[0158]

[0159] Where λ is the decay factor, which can improve the convergence speed by 40% according to actual measurements.

[0160] Introducing a momentum term: Adding a historical gradient weighting term to suppress oscillations.

[0161] ;

[0162] ;

[0163] In the formula, β=0.9 effectively reduces high-frequency oscillations. Let k be the oscillation velocity of the k-th node. Let be the oscillation angle at the k-th node. For historical gradient weighting;

[0164] For example, the calculation of the surface subsidence gradient field;

[0165] 3.1) Spatiotemporal interpolation and gridding;

[0166] Discretize the fused displacement data according to a spatiotemporal grid (recommended resolution: 1km×1km horizontally, 1 day in time), and fill the gaps using Kriging interpolation or inverse distance weighting (IDW).

[0167] 3.2) Gradient field generation;

[0168] Rate calculation: Settlement rate is calculated per grid cell based on time-series displacement data.

[0169]

[0170] In the formula, For the time difference, It is the displacement difference;

[0171] Gradient vector extraction: Calculate the spatial gradient using the Sobel operator or the central difference method.

[0172]

[0173] Characterizes the intensity and direction of spatial variation in settling rate.

[0174] Vector image generation and visualization; GIS platform integration: importing gradient vector fields into QGIS or ArcGIS and overlaying basic geographic information (administrative divisions, topographic elevation); dynamic rendering optimization, including: arrow mapping: vector direction → arrow orientation, gradient magnitude. → Arrow length / color (e.g., red - high settlement, blue - stable zone). Interactive analysis: Supports clicking to query the historical settlement rate curve and confidence interval of any grid point;

[0175] The key implementation tools and verification involved are shown in Table 2.

[0176] Table 2 Key Implementation Tools and Validation

[0177]

[0178] For high-risk areas such as mining areas / landslide areas, it is recommended to cross-validate gradient vector maps with InSAR results to reduce the uncertainty of a single data source; adaptive step size parameters need to be optimized offline using historical data to avoid on-site oscillations and divergences. For example, edge computing and lightweight deployment; on-chip integrated processing unit (MCU+FSM); embedding a finite state machine and Cortex-M3 core within the radar module to achieve in-situ data processing, reducing transmission latency to 59ms; power consumption optimized to 1.7mW, suitable for passive field monitoring scenarios. Cloud-edge collaborative architecture; raw data filtering and feature extraction are completed at the edge, while large-scale gradient field modeling is performed in the cloud, reducing bandwidth usage by 70%; supporting multi-node networking to build regional subsidence monitoring networks, covering an area of ​​square kilometers.

[0179] As can be seen, the core advantages of this invention compared to the prior art are shown in Table 3;

[0180] Table 3. Core advantages of this invention that distinguish it from existing technologies.

[0181]

[0182] Downhole sensing layer: Bionic anchor bolt group sensor, simulating plant root structure, implanted with piezoelectric ceramic and fiber optic grating composite sensing unit, to collect rock stress-acoustic emission signals in real time; for example, the bionic anchor bolt group sensor integration; strain-acoustic emission dual-mode sensing: micro strain gauges (such as fiber optic gratings or vibrating wires) are embedded inside the anchor bolt to monitor stress changes, and piezoelectric ceramic acoustic emission sensors are integrated on the surface of the bolt to capture high-frequency elastic waves (frequency range 50kHz-1MHz) of rock fracture.

[0183] Grid-like layout: Multiple sets of sensors are distributed along the anchor rod axis to form a spatial sensing network, realizing three-dimensional mapping of the stress field (e.g., one measuring point every 20cm). Bionic structural design: Referring to the resonance principle of biological cilia (such as cochlear hair cells), the frequency response characteristics of the sensors are optimized to improve the sensitivity to micro-fracture signals.

[0184] Real-time signal acquisition process;

[0185] graph LR;

[0186] A [rock mass stress / microfracture] --> B (strain sensor → stress electrical signal);

[0187] A --> C (Acoustic emission sensor → elastic wave signal);

[0188] B & C --> D [Signal Preprocessing];

[0189] D --> E [Edge computing node];

[0190] E --> F [Wireless transmission to monitoring platform];

[0191] The real-time acquisition of rock mass stress-acoustic emission signals includes:

[0192] (I) Signal synchronization and noise reduction; Clock synchronization technology: Time synchronization of group sensors is achieved through GPS or high-precision crystal oscillators (error <1μs) to ensure accurate positioning of acoustic emission events. Adaptive filtering: Wavelet transform and EMD decomposition are used to eliminate interference from downhole electromechanical equipment while preserving effective acoustic emission characteristics.

[0193] (II) Edge intelligent processing: Deploy lightweight AI models (such as LSTM networks) on the sensor nodes of the biomimetic anchor group to identify acoustic emission event types (microcrack propagation / rock fracture) in real time and associate them with stress mutation points. Dynamically compressed data: Only upload characteristic parameters of abnormal events (energy, ring count, etc.) to reduce bandwidth requirements.

[0194] (III) Multimodal data fusion; Establishing a stress-acoustic emission correlation model: High-frequency acoustic emission event with sudden stress increase → early warning of local overload of rock mass; Continuous low-energy acoustic emission with gradual stress change → accumulation of creep damage; Among them, data transmission and networking are shown in Table 4.

[0195] Table 4 Data Transmission and Networking

[0196]

[0197] Extreme environment adaptability; Protective design: Silicone sealed housing (IP68) for moisture and explosion protection, pressure resistance range 0-80MPa. Self-powered technology: Piezoelectric energy harvesting device converts anchor vibration into electrical energy, extending battery life. Example: Anti-interference gas monitor: Employs laser spectral compensation algorithm to eliminate the interference of dust... The effect of concentration detection is minimal, with an error rate of <3%.

[0198] For example, hardware-level anti-interference design; dual-wavelength laser scanning technology; principle: two laser beams are emitted simultaneously, one beam is aimed at the gas absorption peak (e.g., CO: 1573nm), 1653nm), another beam scans at a reference wavelength of 317nm without gas absorption. Implementation: Dust scattering interference signals are captured in real time at the reference wavelength. Background noise caused by dust is subtracted through differential calculation, retaining the pure gas absorption signal 426. Case: The DELTER TDLAS ammonia analyzer (BBK-S series) uses this technology to reduce dust interference error to ±0.5%FS17. Optical path structure optimization; long optical path gas chamber design: increasing the number of laser reflections in the gas chamber (e.g., multi-pass cell) to enhance signal strength and offset dust attenuation 326; sheath gas protection system: clean air is introduced at the optical window to form a positive airflow barrier, preventing dust adhesion. Actual measurements show that 0.6MPa pulse purging can reduce dust accumulation at the optical window by 70%. Off-axis integrating cavity enhancement: using high-reflectivity mirrors to construct the optical resonant cavity, the laser is incident at an off-axis angle, dust scattering light is attenuated by the cavity, and the gas absorption signal is amplified by more than 500 times, significantly improving the signal-to-noise ratio.

[0199] For example, algorithm-level dynamic compensation; real-time correction using the Mie scattering model; input parameters: dust concentration (real-time feedback from the built-in particulate sensor), particle size distribution (preset typical values ​​for industrial scenarios). Calculation process:

[0200]

[0201] In the formula, Historical dust scattering concentration This represents the actual measured dust scattering concentration. This is the dust scattering coefficient (related to particle size). For optical path, This refers to the dust concentration.

[0202] Harmonic detection and baseline reconstruction; Second harmonic (2f) extraction: Separating gas absorption signals from low-frequency drift caused by dust through wavelength modulation. Adaptive baseline fitting: Dynamically reconstructing the baseline using signals from the gas-free absorption region to eliminate background drift caused by dust. Multi-physics coupling compensation; Temperature and humidity compensation: The scattering characteristics of dust change after absorbing moisture, requiring simultaneous access to temperature and humidity sensor data and dynamic correction via neural network (e.g., the compensation coefficient increases by 0.15 for every 10% increase in dust humidity). Pressure fluctuation calibration: Dust accumulation may change pipeline air pressure, requiring adjustment of the absorption spectral broadening model in conjunction with pressure sensor data. System-level operation and maintenance strategy; Self-cleaning mechanism; Periodic pulse purging: A 0.6MPa compressed air pulse (lasting 2 seconds) is activated every 5 minutes to remove dust from the optical window. Vibration dust removal: High-frequency micro-vibration of the lens driven by piezoelectric ceramics prevents fine dust deposition. Online diagnosis and calibration; Real-time transmittance monitoring: An alarm is triggered when the transmittance of the optical window is below a threshold (e.g., 80%), prompting cleaning and maintenance. Automatic zero-point calibration: High-purity nitrogen is introduced daily. Correct the baseline to prevent dust accumulation from causing long-term drift.

[0203] Multi-sensor data fusion; linking dust concentration meters (such as laser scattering PM sensors) and gas analyzers to construct a dust-gas interference matrix and improve cross-calibration accuracy. Example: A coal mine inspection robot equipped with a TDLAS methane module and a PM sensor reduced the detection error from ±5% to ±1% after dust compensation. In step S2, the dynamic fracture network evolution model includes: establishing a fractal dimension calculation model for fractures in the goaf.

[0204]

[0205] In the formula, For crack scanning accuracy, To determine the number of fractures, the fracture propagation trend is inverted using microseismic data.

[0206] The coupled seepage-stress equation is:

[0207]

[0208] In the formula, For penetration rate, penetration rate With stress Related, Pore ​​pressure, The saturation level is used to predict the risk of water inrush / gas outburst. For the vector differential operator of the downward gradient, Porosity For a specific moment;

[0209] For example, such as Figure 2 The fracture propagation trend can be inverted using microseismic data, including:

[0210] S201, Data Acquisition and System Deployment: A three-component geophone array is deployed in the target area (such as a mine or fractured well) to form a microseismic monitoring network. The system acquires elastic wave signals (P-wave and S-wave) generated by rock fracturing in real time and records their arrival time, amplitude, and polarity.

[0211] S202, Data preprocessing: Bandpass filtering is applied to the raw seismic records to remove noise interference. Static and dynamic corrections are performed to improve the signal-to-noise ratio and positioning accuracy.

[0212] S203, microseismic event localization; using the P-wave first arrival times from at least four sensors, the source coordinates are calculated using travel-time inversion or grid search methods. Commonly used methods include the Geiger method and improved exhaustive grid search methods, combined with velocity models to improve localization accuracy.

[0213] S204, focal mechanism and parameter inversion; waveform type is determined through polarization analysis, and focal mechanism parameters (strike, dip, slip angle) are extracted; a multi-scan combined waveform inversion algorithm is used to simultaneously optimize spatial location and focal mechanism, and identify shear rupture direction. Microseismic energy, apparent stress, volumetric potential, and other indicators are calculated to reflect rupture intensity and extent.

[0214] For example, polarization analysis is used to determine waveform type, including: polarization analysis is used to distinguish between P-waves (compression waves, linear polarization) and S-waves (shear waves, elliptical / complex polarization). The core formula is based on the eigenvalue decomposition of the covariance matrix.

[0215] 1) Construction of the three-component data matrix: Construct a data matrix from the three components (vertical Z, north-south N, and east-west E) of the data segments recorded by the seismograph.

[0216]

[0217] 2) Calculation of covariance matrix:

[0218]

[0219] In the formula, superscript Let be the transpose matrix, and be the number of sampling points within the time window;

[0220] 3) Eigenvalue decomposition: Solve for the eigenvalues ​​λ1≥λ2≥λ3 and the corresponding eigenvectors v1,v2,v3 of the covariance matrix; where, for P-wave determination, the eigenvalues ​​satisfy λ1≫λ2≈λ3;

[0221] The polarization direction, the eigenvector v1 corresponding to the maximum eigenvalue λ1 indicates the wave propagation direction; among them, for S-wave determination, the eigenvalues ​​satisfy λ1≈λ2>λ3;

[0222] The polarization complexity is high, and combined with the axial ratio of the polarization ellipse, the axial ratio = λ2 / λ1.

[0223] For example, the formula for extracting source mechanism parameters; source mechanism parameters need to be obtained through the P-wave initial motion polarity or full waveform inversion of multiple stations, and polarization analysis is only a preprocessing step.

[0224] (1) P-wave initial motion polarity method (applicable to simple faults); Steps: Identify the direction of P-wave initial motion (upward / compression, downward / expansion). Mark the initial motion symbols of each station in the lower hemisphere projection of the source sphere. Parameter calculation: By finding two orthogonal nodal planes (the fault plane and the auxiliary plane), the initial motion sign is positioned in each of the four quadrants. The parameters are defined geometrically by the nodal planes.

[0225] Strike (φ): The azimuth of the line of intersection of the fault and the horizontal plane relative to true north (range [0°, 360°)).

[0226] Dip angle (δ): The angle between the fault plane and the horizontal plane (range (0°, 90°)).

[0227] Sliding angle (λ): The angle between the sliding direction and the direction of travel (range [-180°, 180°] or [0°, 360°)).

[0228]

[0229] Where λ≈0°, left-lateral strike-slip; λ≈±180°, right-lateral strike-slip; λ≈+90°, reverse fault; Normal fault;

[0230] (2) Waveform inversion method (exact solution, applicable to complex mechanisms); Minimize the objective function:

[0231] Using observed waveforms Compared with theoretical waveforms The residuals are solved by grid search or iterative optimization:

[0232]

[0233] In the formula, For model parameters, ;

[0234] CAP (Cut and Paste) method: The P-wave and surface wave bands are separated and inverted separately using the CAP method, followed by weighted fitting.

[0235]

[0236] In the formula, The weighted fitting values ​​are obtained after inversion of the P-wave and surface wave bands respectively. All are weighting coefficients. To separate the observed waveform residuals from P-wave inversion, To separate the theoretical waveform residuals from P-wave inversion, For the observed waveform residuals of the separated surface wave inversion, The theoretical waveform residual for the separation of surface waves inversion;

[0237] For example, microseismic energy The calculation formula is based on the vibration velocity decay model:

[0238]

[0239] In the formula, The peak vibration velocity at the measuring point is obtained directly through the sensor. The distance from the seismic source to the sensor. This is an energy characteristic coefficient, which is related to the characteristics of the site medium and needs to be calibrated through blasting experiments (e.g., regression analysis of blasting data).

[0240] Application: The higher the energy, the greater the kinetic energy released during rupture, indirectly reflecting the rupture scale.

[0241] Seismic Moment The physical quantity characterizing the fault slip intensity is calculated based on waveform inversion:

[0242]

[0243] In the formula, The density of the medium is (g / cm³). This refers to the seismic wave velocity (P-wave or S-wave, m / s). The low-frequency level of the seismic source spectrum was obtained through spectral analysis. This is the radiation mode coefficient (usually taken as 0.6). This is a directional correction factor;

[0244] Related indicators: seismic moment and moment magnitude The relationship is:

[0245]

[0246] Apparent stress The average stress consumed per unit area of ​​fault slip:

[0247]

[0248] The higher the apparent stress, the stronger the brittle fracture of the rock mass, and it is often used to warn of rockburst risk.

[0249] Body change The expression representing the volumetric change caused by fault slip is:

[0250]

[0251] In the formula, Shear modulus (Pa) is a property of the rock medium (typically 30-70 GPa for granite). Volumetric potential is directly related to fault slip distance and is used to quantify the extent of fracture (e.g., the scale of rock strata displacement in mines).

[0252] Clustering of a large amount of microseismic event point cloud data, the main fracture surface and branch structure are fitted. Combined with prior information such as the geostress field and lithological boundaries, the fracture propagation path and connectivity are deduced. In step S3, a biomimetic intelligent early warning mechanism is implemented; a bee colony optimization neural network (BSO-NN) is constructed: Input layer: 12-dimensional heterogeneous data such as surface deformation, gas concentration, microseismic energy, and fracture density; Hidden layer: a dynamic weight adjustment module, simulating the bee foraging path optimization process, improving response speed by 40%.

[0253] An example is the implementation principle of the dynamic weight adjustment module;

[0254] (1) Mapping of bee foraging behavior; Leading bee (weight exploration): Simulates bees' search for high-quality nectar sources, finding the weight combination that maximizes the accuracy of the risk model in the weight space. Follower bee (weight optimization): Further optimizes the weight allocation locally based on the current optimal weight combination (analogous to nectar source quality). Scout bee (breaking out of local optima): When the weight combination has not been updated for a long time, randomly generates new weights to avoid model stagnation. Adaptive mechanism: Dynamically adjusts the search range through Levy Flight and variable spiral factor to balance global exploration and local development capabilities.

[0255] (2) Correlation between weights and risk indicators; using 12-dimensional heterogeneous data such as surface deformation, gas concentration, microseismic energy, and fracture density as input features, the weight module dynamically allocates the importance of each parameter in the risk model. An example is the implementation process of risk classification and early warning;

[0256] Phase 1: Data acquisition and feature extraction; fusion of multi-source monitoring data, integrating multi-dimensional data such as surface deformation, gas concentration, microseismic energy, and fracture density (12 dimensions in total) to generate time series features; screening of key risk factors, using Pearson correlation coefficient analysis to analyze the correlation between indicators and disasters, and eliminating redundant parameters.

[0257] Phase 2: Construction of the dynamic risk grading model; initial weight allocation; initial weights are generated using Tent chaotic mapping based on historical disaster data to avoid random bias. Iterative optimization of weights (core of the Bee Algorithm);

[0258] # Pseudocode: Weight optimization process;

[0259] for epoch in max_iterations:

[0260] Lead bee updates weights → Calculate fitness (model prediction error);

[0261] Follow the bee to select high fitness weight → neighborhood search optimization;

[0262] If the number of attempts to improve the weight exceeds the limit, the scout bee will randomly reset the weight.

[0263] Save the globally optimal weight combination;

[0264] ``` ```;

[0265] The threshold settings for the different levels are shown in Table 5.

[0266] Phase 3: Real-time Early Warning and Response; Early warning triggering mechanism: when the model outputs a risk level that jumps, a multi-level response is triggered (see Table 6); Dynamic threshold adjustment: weights are recursively updated based on new monitoring data to achieve model self-calibration. It can be seen that the dynamic adaptability of this invention, with weights automatically adjusting according to environmental changes, is superior to fixed-weight models (such as the AHP method). It resists local optima; the reconnaissance bee mechanism effectively avoids the model getting trapped in suboptimal solutions, improving early warning accuracy. It is computationally efficient, with an algorithm complexity of only O(N), supporting minute-level risk updates (compared to hours for traditional finite element simulations).

[0267] Table 5. Tiered Early Warning Strategy

[0268]

[0269] Table 6. Risk Classification Standards for Goaf Area Disasters in Multi-Level Responses

[0270]

[0271] For example, the Bee Swarm Optimized Neural Network (BSO-NN) is a hybrid model that combines swarm intelligence optimization algorithms with neural networks. It mainly optimizes neural network parameters by simulating bee foraging behavior, thereby improving its performance and generalization ability.

[0272] The technical features are as follows:

[0273] (1) Swarm intelligence optimization mechanism; role division: The algorithm simulates three roles in a bee colony: "leader bee" (exploring new solutions), "follower bee" (developing high-quality solutions), and "scout bee" (escaping local optima). Collaborative search: A balance between global and local search is achieved through information sharing (such as "dance zone" transmitting nectar source quality information).

[0274] (2) Neural network parameter optimization; Optimization objective: To address the sensitivity of hyperparameters such as initial weights, thresholds, and learning rates in neural networks, the bee colony algorithm searches for the optimal parameter combination through multiple iterations. Avoiding local optima: A scout bee mechanism is introduced, which randomly resets the search direction when the solution has not been improved for a long time, thereby enhancing the global search capability.

[0275] (3) Hybrid architecture design; Global + Local optimization: The bee colony algorithm is responsible for global parameter search, and the BP neural network is responsible for local gradient fine-tuning, forming a two-layer optimization structure. It includes an input layer and a hidden layer; Adaptive adjustment: The search range (such as the radius of the food source) is dynamically adjusted to adapt to the optimization needs of different training stages.

[0276] (4) High efficiency in parallelism; no backpropagation variant: partially improved model Supports local learning without backpropagation, accelerating the training process. The bee colony optimization phase includes:

[0277] Step 1: Initialization, randomly generate bee colony positions (i.e., combinations of neural network parameters), with each position representing a potential solution.

[0278] Step 2: Leader bees explore, each leader bee generates new solutions in its neighborhood (e.g., adjusting weights), and evaluates the quality of the solutions using a fitness function (e.g., classification accuracy).

[0279] Step 3: Follower bee development. The follower bee focuses on high-quality solutions based on the quality probability of the solution (roulette wheel selection) and further conducts local search.

[0280] Step 4: The scout bee escapes. If a solution fails to improve after multiple attempts, it is abandoned and a new solution is randomly generated to avoid getting trapped in a local optimum.

[0281] The neural network training phase includes:

[0282] Parameter injection: Inject the optimal parameters of bee colony optimization into the neural network (such as initial weights and learning rate).

[0283] Supervised learning: Perform standard forward and backward propagation based on the optimized parameters to minimize the loss function.

[0284] Collaborative working mechanism; dynamic feedback: After each training of the neural network, its output error is fed back to the bee colony algorithm as a fitness value to guide the next round of search.

[0285] Termination condition: Stop when the maximum number of iterations is reached or the quality of the solution is stable, and output the optimal network model.

[0286] As can be seen from the above embodiments, the accuracy is improved: the false alarm rate of surface deformation monitoring is reduced by 67% (compared to patent CN118794484A); the gas anomaly detection time is advanced to more than 72 hours (traditional methods <24 hours).

[0287] Cost optimization: The lifespan of bionic sensors is increased to 5 years (compared to about 2 years for traditional equipment); data transmission adopts a mining-grade Mesh self-organizing network, reducing wired deployment costs by 30%.

[0288] This invention overcomes three major technical shortcomings of existing technologies: it solves the data fusion bottleneck by pioneering a multi-field coupling model of "geology-engineering-fluid"; it addresses weak anti-interference capabilities through biomimetic sensing design and laser spectral compensation; and it resolves static early warning by employing a dynamic fracture network and biomimetic algorithm in a collaborative iteration. In the above embodiments, the descriptions of each embodiment have their own emphasis; for parts not detailed or described in a particular embodiment, please refer to the relevant descriptions in other embodiments.

[0289] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention and within the spirit and principles of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A method for forecasting and early warning of geological disasters in coal mining goaf areas, characterized in that, The method includes the following steps: S1, deploys an aboveground monitoring layer to calculate the surface subsidence gradient field in real time, and deploys a downhole sensing layer to collect rock stress-acoustic emission signals in real time. Gas concentration data; S2. The fractal dimension of the fractures in the goaf is obtained by using a dynamic fracture network evolution model, and the fracture propagation trend is inverted by microseismic data to obtain the relationship between seepage and stress in the goaf. S3 inputs the acquired information into the bee colony optimization neural network BSO-NN, and obtains graded early warning information on geological disasters in the goaf of the coal mining area through a biomimetic intelligent early warning mechanism.

2. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 1, characterized in that, In step S1, the improved millimeter-wave interferometric radar operates at a frequency of 80 GHz in the well monitoring layer. The real-time calculation of the surface subsidence gradient field by integrating BeiDou / GNSS data includes: introducing deep learning YOLOv10 to automatically identify the subsidence funnel region, and combining a two-dimensional Gaussian model to invert nonlinear deformation to solve the phase aliasing of large gradient subsidence. And the calculation of the surface subsidence gradient field.

3. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 2, characterized in that, The construction of a two-dimensional Gaussian model includes: For the inversion of settling funnel parameters, a two-dimensional Gaussian model is constructed for each funnel region identified by YOLOv10, with the objective function of minimizing the angle deviation. The expression is as follows: In the formula, The area of ​​the funnel region. It is horizontal. Vertically, As the central location, To control the range of deformation; The deformation gradient is inverted by fitting the parameters of a two-dimensional Gaussian model through nonlinear optimization.

4. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 2, characterized in that, Calculation of the surface subsidence gradient field includes: Spatiotemporal interpolation and gridding: The fused displacement data is discretized according to a spatiotemporal grid, and Kriging interpolation or inverse distance weighting (IDW) is used to fill the gaps. Gradient field generation includes: calculating settlement rate based on time-series displacement data, per grid cell, with the expression: In the formula, For the time difference, It is the displacement difference; Gradient vector extraction involves calculating the spatial gradient using the Sobel operator or the central difference method, expressed as: Characterizes the intensity and direction of spatial variation in settling rate.

5. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 1, characterized in that, In step S1, the downhole sensing layer includes: a biomimetic anchor bolt group sensor, which simulates the structure of plant roots and is embedded with a piezoelectric ceramic and fiber optic grating composite sensing unit to collect rock stress-acoustic emission signals in real time; The biomimetic anchor bolt group sensor includes: Strain-acoustic emission dual-mode sensing: fiber optic gratings or vibrating wire micro strain gauges are embedded inside the anchor bolt to monitor stress changes, and piezoelectric ceramic acoustic emission sensors are integrated on the surface of the bolt to capture high-frequency elastic waves from rock fracture. The grid-like layout, with multiple sets of sensors distributed along the anchor bolt axis, forms a spatial sensing network to achieve three-dimensional mapping of the stress field. The biomimetic structural design references the resonance principle of biological cilia to design a biomimetic anchor group sensor. Real-time acquisition of rock mass stress-acoustic emission signals includes: (1) Signal synchronization and noise reduction: Time synchronization of group sensors is achieved through GPS or high-precision crystal oscillator; wavelet transform and EMD decomposition are used to eliminate interference from downhole electromechanical equipment while retaining effective acoustic emission characteristics; (2) Edge intelligent processing: deploy a lightweight LSTM network on the sensor nodes of the bionic anchor group to identify acoustic emission event types and associate stress mutation points in real time; dynamically compress data and upload only the characteristic parameters of abnormal events; (3) Multimodal data fusion; establish a stress-acoustic emission correlation model, and issue a warning of local overload of rock mass when stress suddenly increases and high-frequency acoustic emission events are added; when stress gradually changes and continuous low-energy acoustic emission is added, creep damage accumulation is detected. The real-time acquisition of rock mass stress-acoustic emission signals is also corrected in real time using the Mie scattering model, including: Input parameters: dust concentration, particle size distribution; calculation process: In the formula, Historical dust scattering concentration This represents the actual measured dust scattering concentration. The dust scattering coefficient is... For optical path, This refers to the dust concentration.

6. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 1, characterized in that, In step S2, the dynamic fracture network evolution model includes: Establish a model for calculating the fractal dimension of fractures in goaf areas: In the formula, For crack scanning accuracy, To determine the number of fractures, the fracture propagation trend is inverted using microseismic data. The coupled seepage-stress equation is: In the formula, For penetration rate, penetration rate With stress Related, Pore ​​pressure, The saturation level is used to predict the risk of water inrush / gas outburst. For the vector differential operator of the downward gradient, Porosity For a specific moment; Inversion of fracture propagation trends using microseismic data includes: S201, Data Acquisition and System Deployment: A three-component geophone array is deployed in the target area to form a microseismic monitoring network; the system acquires P-wave and S-wave elastic wave signals generated by rock fracture in real time, recording the time, amplitude and polarity. S202, Data preprocessing; bandpass filtering of the raw seismic record; static and dynamic corrections; S203, Microseismic event location; using the P-wave first arrival times of at least 4 sensors, the source coordinates are calculated using travel-time inversion or grid search methods; S204, focal mechanism and parameter inversion; waveform type is determined by polarization analysis, and focal mechanism parameters, including strike, dip, and slip angle, are extracted; a multi-scan combined waveform inversion algorithm is used to simultaneously optimize spatial location and focal mechanism, and identify shear rupture direction; microseismic energy, apparent stress, and volumetric potential index are calculated to reflect rupture intensity and extent.

7. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 6, characterized in that, Polarization analysis determines waveform type, including: Three-component data matrix construction: A data matrix is ​​constructed from the three component data segments recorded by the seismograph: vertical Z, north-south N, and east-west E. Covariance matrix calculation: In the formula, superscript Let be the transpose matrix, and be the number of sampling points within the time window; Eigenvalue decomposition: Solve for the eigenvalues ​​λ1≥λ2≥λ3 and the corresponding eigenvectors v1,v2,v3 of the covariance matrix; where, for P-wave determination, the eigenvalues ​​satisfy λ1≫λ2≈λ3; The polarization direction, the eigenvector v1 corresponding to the maximum eigenvalue λ1 indicates the wave propagation direction; among them, for S-wave determination, the eigenvalues ​​satisfy λ1≈λ2>λ3; The polarization complexity is high, and combined with the axial ratio of the polarization ellipse, the axial ratio = λ2 / λ1.

8. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 7, characterized in that, The focal mechanism parameters were obtained by inverting the P-wave initial polarity or full waveform from multiple stations; The P-wave initial polarity method includes: Identify the P-wave initial motion direction; mark the initial motion symbols of each station in the lower hemisphere projection of the source sphere; calculate parameters by finding two orthogonal nodal planes so that the initial motion symbols are located in the four quadrants; the parameters are defined geometrically by the nodal planes. The strike φ, the azimuth of the line of intersection of the fault and the horizontal plane relative to true north, ranging from [0°, 360°); Dip angle δ, the angle between the fault plane and the horizontal plane, ranging from 0° to 90°; Sliding angle λ is the angle between the sliding direction and the direction of travel, ranging from [-180°, 180°] to [0°, 360°]. Where λ≈0°, left-lateral strike-slip; λ≈±180°, right-lateral strike-slip; λ≈+90°, reverse fault; Normal fault; The full waveform inversion method includes: Minimize the objective function using observed waveforms Compared with theoretical waveforms The residuals are solved by grid search or iterative optimization: In the formula, For model parameters, ; P-wave and surface wave bands were separated and inverted separately using the CAP method, followed by weighted fitting: In the formula, The weighted fitting values ​​are obtained after inversion of the P-wave and surface wave bands respectively. All are weighting coefficients. To separate the observed waveform residuals from P-wave inversion, To separate the theoretical waveform residuals from P-wave inversion, For the observed waveform residuals of the separated surface wave inversion, The theoretical waveform residual for the separation of surface waves inversion; Micro-vibration energy The calculation formula is based on the vibration velocity decay model: In the formula, The peak vibration velocity at the measuring point is obtained directly through the sensor. This represents the distance from the seismic source to the sensor. This is an energy characteristic coefficient, which is related to the characteristics of the site medium. Seismic Moment The physical quantity characterizing the fault slip intensity is calculated based on waveform inversion: In the formula, For the density of the medium, For seismic wave velocity, The low-frequency level of the seismic source spectrum was obtained through spectral analysis. For radiation mode coefficients, This is a directional correction factor; Related indicators: seismic moment and moment magnitude The relationship is: Apparent stress The average stress consumed per unit area of ​​fault slip: Body change The expression representing the volumetric change caused by fault slip is: In the formula, This is the shear modulus.

9. The method for forecasting and early warning of geological disasters in coal mining goaf areas according to claim 1, characterized in that, In step S3, the bee colony optimization neural network BSO-NN includes: Input layer: multidimensional heterogeneous data on surface deformation, gas concentration, microseismic energy, and fracture density; Hidden layer, dynamic weight adjustment module, simulates the process of bee foraging path optimization; The implementation process of risk classification and early warning includes: Phase 1: Data acquisition and feature extraction; fusion of multi-source monitoring data, integrating multi-dimensional heterogeneous data such as surface deformation, gas concentration, microseismic energy, and fracture density to generate time series features; screening of key risk factors, using Pearson correlation coefficient analysis to analyze the correlation between indicators and disasters, and eliminating redundant parameters; Phase 2: Construction of dynamic risk classification model; initial weight allocation; based on historical disaster data, generate initial weights using Tent chaotic mapping and set classification thresholds; Phase 3: Real-time early warning and response; early warning triggering mechanism, which triggers multi-level response when the model outputs a risk level that jumps; dynamic threshold adjustment, which recursively updates weights based on new monitoring data to achieve model self-calibration.

10. A geological disaster forecasting and early warning system for goaf areas in coal mines, characterized in that, This system implements the geological disaster forecasting and early warning method for goaf areas in coal mining areas as described in any one of claims 1-9. The system includes: The data acquisition module is used to deploy an aboveground monitoring layer to calculate the surface subsidence gradient field in real time; and to deploy a downhole sensing layer to acquire rock mass stress-acoustic emission signals in real time, as well as to acquire... Gas concentration data; The goaf seepage-stress relationship acquisition module is used to obtain the fractal dimension of the goaf fracture using a dynamic fracture network evolution model, invert the fracture propagation trend through microseismic data, and obtain the goaf seepage-stress relationship. The graded early warning module is used to input the acquired information into the bee colony optimization neural network BSO-NN, and obtain graded early warning information on geological disasters in the goaf of coal mining areas through a biomimetic intelligent early warning mechanism.

Images