Coal mine geological modeling system based on geophysical prospecting data processing and machine learning

Abstract

The invention discloses a coal mine geological modeling system based on geophysical prospecting data processing and machine learning, and relates to the technical field of coal mining. Comprising a multi-modal sensing module for fusing multi-source heterogeneous data and determining a fusion feature value; the physical modeling module is used for obtaining a pressure residual error corresponding to a wave equation regularization layer by fusing characteristic values, adjusting the trend of the established combined learning model, and obtaining a final combined learning model; and the safety decision module obtains a risk probability through the final joint learning model and real-time data, and triggers an early warning signal according to the risk probability and an early warning risk threshold. The obtained multi-modal data are fused, the limitation of a traditional single data source is overcome, meanwhile, data with different sampling rates and spatial-temporal scales are effectively aligned through the dynamic time warping and space mapping technology, and the real-time performance and accuracy of data fusion are further improved.

Description

Technical Field

[0001] The present invention relates to the technical field of coal mining, and in particular to a coal mine geological modeling system based on geophysical data processing and machine learning. Background Art

[0002] The mining depth of my country's coal resources is extending at a rate of 10-15 meters per year. Complex structures such as faults and collapse columns are common in deep coal seams. The cost of a single hole based on traditional drilling sampling exceeds 800,000 yuan per 100 meters, and the drilling spacing is usually greater than 200 meters, which makes it difficult to meet the requirements of the "Coal Mine Safety Regulations" for the detection accuracy of hidden disaster-causing bodies. Although geophysical exploration technology (such as 3D seismic exploration and transient electromagnetic method) can achieve large-scale detection, its data interpretation is highly dependent on the experience of engineers, especially in areas with significant anisotropy of coal-bearing strata, which will increase the uncertainty of the inversion results, resulting in insufficient modeling accuracy.

[0003] The Chinese invention patent with the publication number CN105866836A discloses a method for predicting the full-process geological exploration of a mine using three-dimensional seismic exploration methods, which includes the following steps: using a three-dimensional seismic exploration method to delineate whether there are concealed geological structures in the mineable coal layer; using the actual tunnel exposure and drilling results as constraints, the results of the three-dimensional seismic exploration are verified and analyzed and evaluated in real time; the verified and analyzed geological information of the mineable coal layer is integrated into the three-dimensional seismic exploration results, and the exploration area results that do not meet the detection accuracy requirements are reprocessed, reinterpreted and reapplied using an iterative feedback method. The constraints are constantly changing, and finally the results of the full-process geological analysis of the three-dimensional seismic in the coal mine are formed. This method can effectively improve the accuracy of three-dimensional seismic processing and interpretation, and can predict and identify the geological structure within the scope of mining influence to the maximum extent, which is convenient for the layout of the tunnel mining area, and can scientifically and effectively avoid and reduce production safety accidents.

[0004] The above-mentioned and similar geological exploration methods mostly rely on static historical data to construct the initial three-dimensional model during the geological survey process. It is difficult to effectively integrate the multi-source heterogeneous geophysical data obtained in real time during the mining process (such as microseismic monitoring, logging while drilling, electromagnetic wave CT, etc.). At the same time, due to the characteristics of high noise, multiple sampling rates and mismatch of time and space scales of such data, the model update cycle will be as long as several hours to several days, resulting in a significant time and space deviation from the advancement speed of several meters per hour of the mining face. Summary of the invention

[0005] The purpose of the present invention is to provide a coal mine geological modeling system based on geophysical data processing and machine learning to solve the problems raised in the above background technology.

[0006] To achieve the above object, the present invention provides the following technical solutions: A coal mine geological modeling system based on geophysical data processing and machine learning, including:

[0007] A multi-modal perception module, which obtains the difference magnitude corresponding to each piece of the multi-source heterogeneous data through the collected multi-source heterogeneous data, and fuses the multi-source heterogeneous data through the difference magnitude to determine a fusion feature value;

[0008] SA1: Obtain multi-source heterogeneous data streams: Obtain multi-source heterogeneous data within the advancing direction of the working face through microseismic sensors, drill pipe devices, and electromagnetic wave CT pairs;

[0009] SA2: Perform dynamic time warping: Obtain the coordinates and time corresponding to the effective energy through a microseismic sensor, obtain the interface position and measurement time of the well logging through a drill pipe device, and align the coordinates corresponding to the effective energy and the interface position, and the time corresponding to the effective energy and the measurement time;

[0010] Step SA3: Obtain a three-dimensional fusion feature volume: Obtain corresponding interpolation data through the effective energy obtained by a microseismic sensor, the gamma ray intensity while drilling, and the CT dielectric constant of the electromagnetic wave, and obtain a fused normalized value according to the interpolation data;

[0011] A physical modeling module, which obtains the pressure residual corresponding to the wave equation regularization layer through the fusion feature value, and adjusts the trend of the established joint learning model according to the pressure residual to obtain a final joint learning model;

[0012] A safety decision-making module, which obtains the risk probability corresponding to each piece of the real-time data through the final joint learning model and the real-time data, and triggers an early warning signal according to the risk probability and the early warning risk threshold.

[0013] Furthermore, obtaining the multi-source heterogeneous data within the advancing direction of the working face includes:

[0014] SA1.1: Deploy microseismic monitoring nodes: Obtain effective energy signals through microseismic sensors and a preset energy threshold, specifically:

[0015]

[0016] Where: E is the signal energy of the microseismic event, t1 is the start point of the integration time window, t2 is the end point of the integration time window, a x is the instantaneous acceleration value of the microseismic sensor in the x orthogonal direction, a y is the instantaneous acceleration value of the microseismic sensor in the y orthogonal direction, a z is the instantaneous acceleration value of the microseismic sensor in the z orthogonal direction;

[0017] SA1.2: Set the logging-while-drilling unit: Correct the azimuth deviation of the detector through the attitude angle of the drill bit and the coordinate transformation matrix, specifically:

[0018]

[0019] Where: x′ is the roadway driving direction, y′ is the lateral direction perpendicular to the roadway, z′ is the gravity direction perpendicular to the downward direction, and R z (ψ) is the yaw angle rotation matrix about the Z axis, ψ is the yaw angle of the drill bit, and R y (θ) is the pitch angle rotation matrix about the y axis, θ is the pitch angle of the drill bit, x is the drilling direction along the drill pipe axis, y is the lateral direction perpendicular to the drill pipe axis, and Z is the vertical direction perpendicular to the drill pipe axis;

[0020] SA1.3: Configure the electromagnetic wave CT pair: Obtain the true coal and rock dielectric constant through the electromagnetic wave transmitter and the array receiver, specifically:

[0021]

[0022] Where: ε r is the true coal and rock dielectric constant, E m is the measured electromagnetic field strength, E s is the scattered field strength caused by the metal support structure, E t is the theoretical field strength without metal interference, and ||·|| 2 is the square of the norm.

[0023] Furthermore, align the coordinates and interface positions corresponding to the effective energy, and the time corresponding to the effective energy and the measurement time, including:

[0024] SA2.1: Data timestamp alignment: Obtain the total cost function of time according to the time series of the effective energy signal and the time series of the logging, specifically:

[0025] Tot(i, j) = Tim(i, j) + λ·IMU(i, j)

[0026] Where: Tot(i, j) is the total alignment cost, Tim(i, j) is the time difference cost, λ is the weight coefficient, IMU(i, j) is the attitude difference cost, i is the time sequence index of the microseismic event, and j is the time sequence index of the logging-while-drilling;

[0027] SA2.2: Spatial mapping: Obtain the three-dimensional coordinates in the set roadway coordinate system through the UWB tag and the mobile tag, and adjust the position of each USB base station according to the three-dimensional coordinates.

[0028] Further, the positions of each USB base station are adjusted, including:

[0029] SA2.2.1: Set up the roadway coordinate system: Through the UWB tag, take the position of the cutting head of the roadheader as the origin of the roadway coordinate system, take the advancing direction of the roadheader as the X-axis, the vertically upward direction of the roadheader as the Y-axis, and the horizontal transverse direction of the roadheader as the Z-axis to establish the roadway coordinate system;

[0030] SA2.2.2: Correct the position of the UWB base station: UWB base stations are arranged at intervals in the roadway, mobile tags are set in the microseismic sensors and the logging-while-drilling equipment. At the same time, in the roadway coordinate system, obtain the three-dimensional coordinates of each mobile tag, and obtain the Euclidean distance between the three-dimensional coordinates and the UWB base station. And adjust the position of the UWB base station according to the magnitude relationship between the Euclidean distance and the preset error. Specifically:

[0031] When the Euclidean distance is not less than the preset error, adjust the position of the UWB base station until the Euclidean distance is less than the preset error. Otherwise, the position of the UWB base station is not adjusted.

[0032] Further, obtain the fused normalized value, including:

[0033] SA3.1: Obtain Kriging interpolation: Through the three-dimensional coordinates of each data point, perform interval grouping, and determine the average semivariogram value corresponding to each interval grouping according to the attribute differences corresponding to all data points in each interval grouping;

[0034] SA3.2: Obtain interpolation weights: According to the position of the center point of the target voxel and the attribute values of all data points within the preset range, obtain the weight coefficient of each data point, and determine the attribute interpolation result of the center point of the target voxel according to the weight coefficient;

[0035] SA3.3: Multimodal feature fusion: Perform normalization processing on the attribute interpolation result, and obtain the fusion feature value according to the normalized attribute interpolation result. Specifically:

[0036] L = α1·E nor +α2·ζ nor +α3·ε nor

[0037] Where: L is the fusion feature value, α1 is the weight coefficient of the effective energy signal, E nor is the normalized value of the effective energy signal, α2 is the weight coefficient of the gamma ray, ζ nor is the normalized value of the gamma ray, α3 is the weight coefficient of the coal-rock dielectric constant, ε noris the normalized value of the dielectric constant of coal and rock.

[0038] Further, determining the average semivariogram value corresponding to each of the interval groups includes:

[0039] SA3.1.1: Conduct semivariogram modeling: According to the three-dimensional coordinates of each data point in the roadway coordinate system, obtain the three-dimensional Euclidean distance between the data points, specifically:

[0040]

[0041] where: h is the straight-line distance between two data points, x m is the X-axis coordinate of the mth data point in the roadway coordinate system, x n is the X-axis coordinate of the nth data point in the roadway coordinate system, y m is the Y-axis coordinate of the mth data point in the roadway coordinate system, y n is the Y-axis coordinate of the nth data point in the roadway coordinate system, z m is the Z-axis coordinate of the mth data point in the roadway coordinate system, z n is the Z-axis coordinate of the nth data point in the roadway coordinate system;

[0042] SA3.1.2: Conduct distance grouping: According to the three-dimensional Euclidean distance between the data points and a preset distance, conduct interval grouping;

[0043] SA3.1.3: Obtain the average semivariogram value: According to the attribute difference between the data points in each of the interval groups, obtain the average semivariogram value corresponding to each interval group, specifically:

[0044]

[0045] where: γ(h) is the semivariogram value, N(h) is the number of pairs of data points near the lag distance h, Z(X r ) is the attribute value corresponding to the position X r , Z(X r +h) is the attribute value corresponding to another point position at a distance h from the position X r , X r is the three-dimensional coordinate of the rth data point in the interval group in the roadway coordinate system, r is the index of the data point in the interval group, and h is the straight-line distance between two data points.

[0046] Further, determining the attribute interpolation result of the center point of the target voxel includes:

[0047] SA3.2.1: Construct the Kriging equations: Based on the center point position of the target voxel and all data points within the preset range, establish the Kriging equations through the semivariogram model to determine the weight coefficient corresponding to each data point. The specific Kriging equations are as follows:

[0048]

[0049] where: ω d is the d-th weight coefficient, γ(h dp ) is the semivariogram value between the d-th neighboring data point and the p-th neighboring data point, μ is the Lagrange multiplier, γ(h d0 ) is the semivariogram value between the d-th neighboring data point and the center point of the target voxel, d is the index of the weight coefficient and the neighboring data point, h dp is the straight-line distance between the d-th neighboring data point and the p-th neighboring data point, h d0 is the straight-line distance between the d-th neighboring data point and the center point of the target voxel, and n is the number of neighboring data points participating in the interpolation;

[0050] SA3.2.2: Obtain the attribute interpolation result: Based on the weight coefficient and attribute value corresponding to the data point, determine the attribute interpolation result of the center point of the target voxel, specifically:

[0051]

[0052] where: R is the attribute interpolation result of the center point of the target voxel, n is the number of neighboring data points participating in the interpolation, d is the index of the weight coefficient and the neighboring data point, ω d is the d-th weight coefficient, and Z d is the attribute value of the d-th neighboring data point.

[0053] Furthermore, obtaining the final joint learning model includes:

[0054] SB1: Obtain the residual: Through the fused eigenvalue, obtain the predicted wave velocity, and based on the predicted wave velocity and the predicted pressure obtained by the neural network model, obtain the wave equation residual, specifically:

[0055]

[0056] where: R w is the wave equation residual, p pred is the predicted pressure field corresponding to the neural network model, T is the discretized time, v pred is the predicted wave velocity, e, f, g are the three-dimensional spatial grid indices, h is the index of the discretized time, is the Laplace operator;

[0057] SB2: Obtaining the model trend: According to the residual of the wave equation, obtain the trend of the joint learning model, and adjust the joint learning model according to the trend. The specific formula for obtaining the trend is as follows:

[0058]

[0059] Where: R w is the residual of the wave equation, p pred is the predicted pressure field corresponding to the neural network model, T is the discretized time, v pred is the predicted wave velocity, is the Laplace operator, and θ is the model trend.

[0060] Furthermore, the specific formula for obtaining the predicted wave velocity is as follows:

[0061]

[0062] Where: v pred is the predicted wave velocity, μ is the Lagrange multiplier, and L is the fusion eigenvalue.

[0063] Furthermore, triggering the warning signal includes:

[0064] SC1: Obtaining the risk threshold: Through the microseismic energy, gas concentration, stress field concentration, and the final joint learning model, obtain the risk probability, specifically:

[0065]

[0066] Where: P risk is the risk probability, β1 is the weight coefficient of the signal energy, β2 is the weight coefficient of the gas concentration, β3 is the weight coefficient of the stress gradient, b is the bias term, E is the signal energy of the microseismic event, C is the gas concentration, is the stress gradient;

[0067] SC2: Triggering the warning signal: Compare the risk probability with the preset risk threshold, and trigger the warning signal according to the comparison result, specifically:

[0068] When the risk probability is greater than the first-level risk threshold but less than the second-level risk threshold, trigger the first-level warning signal. When the risk probability is not less than the second-level risk threshold, trigger the second-level warning signal. Otherwise, do not trigger the warning signal.

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

[0070] First: The present invention performs fusion processing on multi-modal data obtained by microseismic sensors, logging while drilling, and electromagnetic wave CT, overcoming the limitations of traditional single data sources. At the same time, through dynamic time warping and spatial mapping technologies, data with different sampling rates and spatio-temporal scales are effectively aligned, further improving the real-time performance and accuracy of data fusion;

[0071] Second: The present invention is based on semi-variogram modeling of Kriging interpolation, combined with weighted fusion of multi-modal features, to generate a three-dimensional fusion feature volume, solving the interpolation error problem under sparse data and improving the spatial resolution of geological modeling;

[0072] Third: The present invention dynamically corrects the base station position through the Euclidean distance feedback between the UWB base station and the mobile tag, reducing the positioning error caused by roadway environment interference and ensuring the spatial consistency of multi-source data in a unified coordinate system. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a schematic flow chart of the coal mine geological modeling system of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0074] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0075] Refer to Figure 1 , this embodiment provides a coal mine geological modeling system based on geophysical data processing and machine learning. The coal mine geological modeling system includes a multi-modal perception module, a physical modeling module, and a safety decision-making module. Among them, the multi-modal perception module obtains the difference size corresponding to each multi-source heterogeneous data through the collected multi-source heterogeneous data, and fuses the multi-source heterogeneous data through the difference size to determine the fusion feature value. The physical modeling module obtains the wave equation regularization layer according to the fusion feature value, and adjusts the trend of the established joint learning model to obtain the final joint learning model. The safety decision-making module obtains the risk probability corresponding to each real-time data through the final joint learning model and real-time data, and triggers an early warning signal according to the risk probability and the early warning risk threshold.

[0076] In this embodiment, the multi-modal perception module obtains multi-source heterogeneous data in the advancing direction of the mining face through microseismic sensors, drill pipe devices, and electromagnetic wave CT. At the same time, corresponding preprocessing is performed on the obtained multi-source heterogeneous data, and the preprocessed multi-source heterogeneous data is fused to obtain the fusion feature value. Specifically as follows:

[0077] Step SA1: Obtain multi-source heterogeneous data streams. That is, monitor the energy threshold of nodes through microseismic sensors to capture effective energy in real time and filter out the impact noise of picks. At the same time, measure the attitude angle of the drill bit through the gyroscope and accelerometer on the drill pipe device, and correct the azimuth deviation of the detector. At the same time, obtain the measured electromagnetic field strength and the scattered field strength caused by the metal support structure through the electromagnetic wave CT transmitter / receiver pair to correct the dielectric constant of coal and rock. That is to say, through the microseismic sensors, drill pipe devices and electromagnetic wave CT pairs, multi-source heterogeneous data streams within 20 meters in front of the advancing direction of the mining face are collected in real time, specifically as follows:

[0078] Step SA1.1: Deploy microseismic monitoring nodes. That is, install microseismic sensors at intervals on the roof of the mining roadway, and capture effective energy signals in real time through the microseismic sensors and the preset energy threshold, specifically:

[0079]

[0080] Where: E is the signal energy of the microseismic event, t1 is the start point of the integration time window, t2 is the end point of the integration time window, a x is the instantaneous acceleration value of the microseismic sensor in the x orthogonal direction, a y is the instantaneous acceleration value of the microseismic sensor in the y orthogonal direction, a z is the instantaneous acceleration value of the microseismic sensor in the z orthogonal direction.

[0081] In the process of specific implementation, four-component microseismic sensors are installed at intervals of 10 meters along the roof of the mining roadway. The frequency band range of the four-component microseismic sensor is 50 - 500 Hz, and the sampling rate ≥ 1 kHz. Further, the preset energy threshold is set to 10 J, that is, when the signal energy of the microseismic event is not less than the preset energy threshold (i.e., 10 J), the four-component microseismic sensor performs signal acquisition.

[0082] Step SA1.2: Set up the logging-while-drilling unit. That is, set a detector and a resistivity probe at the front end of the drill pipe, and at the same time, measure the attitude angle of the drill pipe end (i.e., the drill bit) in real time through the gyroscope and accelerometer built into the drill pipe. The attitude angle of the drill bit includes the pitch angle and yaw angle of the drill bit. At the same time, according to the obtained pitch angle and yaw angle, and through the coordinate transformation matrix, correct the azimuth deviation of the detector, specifically:

[0083]

[0084] Where: x′ is the roadway advancing direction, y′ is the lateral direction perpendicular to the roadway, z′ is the gravity direction vertically downward, R z$(\psi)$ is the yaw angle rotation matrix about the Z-axis, $\psi$ is the yaw angle of the drill bit, $R$ y $(\theta)$ is the pitch angle rotation matrix about the y-axis, $\theta$ is the pitch angle of the drill bit, $X$ is the drilling direction along the axial direction of the drill pipe, $y$ is the lateral direction perpendicular to the axial direction of the drill pipe, and $Z$ is the vertical direction perpendicular to the axial direction of the drill pipe.

[0085] Step SA1.3: Configure the electromagnetic wave CT pair. Among them, the electromagnetic wave transmitter points vertically to the coal seam and emits electromagnetic waves outward. The array receivers are arranged at intervals to receive electromagnetic waves and obtain the amplitude and phase of the electromagnetic waves.

[0086] Furthermore, based on the roadway support design drawings, a 3D model of the metal bolt / support is established, and the scattering field intensity caused by the metal support structure is obtained through the boundary element method algorithm. Specifically:

[0087]

[0088] Where: is the Laplace operator, $E$ s is the scattering field intensity caused by the metal support structure, and $k$ is the wave number.

[0089] Furthermore, according to the obtained scattering field intensity caused by the metal support structure and the electromagnetic field intensity measured in real time, the true coal and rock dielectric constant is obtained. Specifically:

[0090]

[0091] Where: $\varepsilon$ r is the true coal and rock dielectric constant, $E$ m is the measured electromagnetic field intensity, $E$ s is the scattering field intensity caused by the metal support structure, $E$ t is the theoretical field strength without metal interference, $\|\cdot\|$ 2 is the square of the norm.

[0092] Step SA2: Perform dynamic time warping. That is, capture the source coordinates and time of the effective energy and the interface position and measurement time of the logging with the drill pipe device through the microseismic sensor, and align the two. Specifically as follows:

[0093] Step SA2.1: Align the data timestamps. That is, according to the time series of the effective energy signal captured by the microseismic sensor and the time series of the logging with the drill pipe device, obtain the total cost function between them. Specifically:

[0094] $Tot(i,j)=Tim(i,j)+\lambda\cdot IMU(i,j)$

[0095] Where: Tot(i, j) is the total alignment cost, Tim(i, j) is the time difference cost, λ is the weight coefficient, IMU(i, j) is the attitude difference cost, i is the time series index of the microseismic event, and j is the time series index of the logging-while-drilling.

[0096] Furthermore, by capturing the time series of the effective energy signal with a microseismic sensor and the time series of logging with a downhole tool, the time difference cost between the two is obtained, specifically:

[0097] Tim(i, j) = |t i - t' j |

[0098] Where: Tim(i, j) is the time difference cost, t i is the acquisition time of the i-th microseismic event, t' j is the acquisition time of the j-th logging-while-drilling, i is the time series index of the microseismic event, and j is the time series index of the logging-while-drilling.

[0099] During the specific implementation process, the acquisition time of the i-th microseismic event is 10.5 s, and the acquisition time of the j-th logging-while-drilling is 10.6 s. Therefore, the time difference cost between the two is 0.1 s.

[0100] Furthermore, according to the pitch angle and yaw angle corresponding to the acquisition time t i of the i-th microseismic event of the drill bit, and at the same time combining it with the attitude difference of the drill bit at the acquisition time t' j of the j-th logging-while-drilling, the attitude difference cost between the two is obtained, specifically:

[0101] IMU(i, j) = |θ i - θ j | + |ψ i - ψ j |

[0102] Where: IMU(i, j) is the attitude difference cost, i is the time series index of the microseismic event, j is the time series index of the logging-while-drilling, θ i is the pitch angle of the drill bit obtained at time t i , θ j is the pitch angle of the drill bit obtained at time t' j , ψ i is the yaw angle of the drill bit obtained at time t i , ψ j is the yaw angle of the drill bit obtained at time t' j , t i is the acquisition time of the i-th microseismic event, and t' j is the acquisition time of the j-th logging-while-drilling.

[0103] During the specific implementation process, at time t i the pitch angle of the drill bit obtained is 5°, and at time t' j the pitch angle of the drill bit obtained is 6°. At time t i the yaw angle of the drill bit obtained is 2°, and at time t' j the yaw angle of the drill bit obtained is 3°. Therefore, the attitude difference cost between the two is 2°.

[0104] Furthermore, based on the time difference cost (0.1 s) and attitude difference cost (2°) obtained between the two above, the total cost function between them, that is, the normalized value is 1.1.

[0105] Step SA2.2: Spatial mapping. That is, through the wireless communication between the UWB tag and the mobile tag, obtain the three-dimensional coordinates in the set roadway coordinate system, and according to the obtained three-dimensional coordinates, adjust the position of each USB base station in real time to ensure smooth wireless communication between the UWB tag and the mobile tag. Specifically as follows:

[0106] Step SA2.2.1: Set the roadway coordinate system. That is, install a UWB tag at the position of the cutting head of the roadheader, and use the position of this cutting head of the roadheader as the origin of the roadway coordinate system. It should be noted that as the position of the roadheader moves, whenever the roadheader moves a preset distance (which can be specifically set according to actual needs), the origin of the roadway coordinate system is re-calibrated.

[0107] Furthermore, at the origin position of the roadway coordinate system, take the advancing direction of the roadheader as the X-axis, that is, the X-axis of the roadway coordinate system is along the roadway driving direction, and it is calibrated with a laser pointer or an inertial navigation system. Furthermore, take the vertically upward direction of the roadheader as the Y-axis, that is, the Y-axis of the roadway coordinate system is the vertical direction determined by a gravity sensor and is orthogonal to the X-axis. Furthermore, take the horizontal transverse direction of the roadheader as the Z-axis, that is, the Z-axis of the roadway coordinate system points to the roadway sidewall and is perpendicular to the X-Y plane.

[0108] Step SA2.2.2: Correct the position of the UWB base station. That is, set UWB base stations at intervals in the roadway, and obtain the position coordinates of each UWB base station through a total station. At the same time, set mobile tags in each microseismic sensor and downhole drilling equipment, and through the wireless communication between the mobile tag and the UWB base station, obtain the three-dimensional coordinates of each mobile tag in the roadway coordinate system established in Step SA2.2.1.

[0109] Furthermore, according to the three-dimensional coordinates of each mobile tag and the position coordinates of each UWB base station, the distance between them is obtained, and the distance between them is compared with a preset error. According to the comparison result, the position of the UWB base station is adjusted in real time to ensure smooth wireless communication between the UWB base station and the mobile tag. Specifically:

[0110] When the Euclidean distance between the three-dimensional coordinates of the mobile tag and the position coordinates of the UWB base station is not less than the preset error, the position of the UWB base station with a distance less than the preset error is adjusted until the Euclidean distance between the three-dimensional coordinates of the mobile tag and the position coordinates of the UWB base station is less than the preset error. Otherwise, the position of the UWB base station is not adjusted.

[0111] In this embodiment, the Euclidean distance between the three-dimensional coordinates of the mobile tag and the position coordinates of the UWB base station is specifically:

[0112]

[0113] where: η U is the Euclidean distance between the three-dimensional coordinates of the mobile tag and the position coordinates of the UWB base station, x U is the X-axis coordinate of the mobile tag in the roadway coordinate system, x r is the X-axis coordinate of the UWB base station, y U is the Y-axis coordinate of the mobile tag in the roadway coordinate system, y r is the Y-axis coordinate of the UWB base station, z U is the Z-axis coordinate of the mobile tag in the roadway coordinate system, z r is the Z-axis coordinate of the UWB base station.

[0114] Step SA3: Obtain the three-dimensional fusion feature body. That is, through the discrete points of the microseismic energy field, the gamma ray intensity while drilling, and the CT dielectric constant of the electromagnetic wave, the effective energy signal, gamma ray, and coal-rock dielectric constant corresponding to each data point are interpolated to obtain the interpolation result corresponding to each data point, and according to the interpolation result, the fused normalized value is obtained. Specifically as follows:

[0115] Step SA3.1: Obtain the Kriging interpolation. That is, through the three-dimensional coordinates of each data point, the three-dimensional Euclidean distance between different data points is obtained for interval grouping, and at the same time, according to the attribute difference between the data points in each interval grouping, the average semivariogram value corresponding to each interval grouping is determined. Specifically as follows:

[0116] Step SA3.1.1: Conduct semivariogram modeling. That is, according to the established roadway coordinate system, the three-dimensional coordinates corresponding to each data point are obtained, and through the three-dimensional coordinates corresponding to each data point, the corresponding three-dimensional Euclidean distance is obtained. Specifically:

[0117]

[0118] Where: h is the straight-line distance between two data points, and x m is the X-axis coordinate of the m-th data point in the roadway coordinate system, and x n is the X-axis coordinate of the n-th data point in the roadway coordinate system, and y m is the Y-axis coordinate of the m-th data point in the roadway coordinate system, and y n is the Y-axis coordinate of the n-th data point in the roadway coordinate system, and z m is the Z-axis coordinate of the m-th data point in the roadway coordinate system, and z n is the Z-axis coordinate of the n-th data point in the roadway coordinate system.

[0119] In the process of specific implementation, the position coordinates corresponding to one data point are (5.0 m, 3.0 m, 2.0 m), and the position coordinates corresponding to another data point are (5.5 m, 3.2 m, 2.1 m). Therefore, the straight-line distance between these two data points is 0.55 m.

[0120] Step SA3.1.2: Perform distance grouping. That is, according to the straight-line distance between the two data points obtained in step SA3.1.1, compare it with the preset distance value to group different straight-line distances into intervals.

[0121] In the process of specific implementation, the preset distance is set to 0.5 m. Therefore, in the process of interval grouping, it can be divided into: [0, 0.5 m), [0.5 m, 1.0 m), and [1.0 m, 1.5 m), and so on.

[0122] Step SA3.1.3: Obtain the average semivariogram value. That is, according to the interval grouping divided in step SA3.1.2, obtain the attribute difference between data points in each interval grouping, and through the attribute difference between data points, obtain the average semivariogram value corresponding to each interval grouping. Specifically:

[0123]

[0124] Where: γ(h) is the semivariogram value, N(h) is the number of pairs of data points near the lag distance h, Z(X r ) is the attribute value corresponding to the position X r at, Z(X r +h) is the attribute value corresponding to the position of another point at a distance h from the position X r , and X r$P_r$ is the three-dimensional coordinate of the $r$-th data point in the roadway coordinate system for interval grouping, $r$ is the index of the data point in the interval grouping, and $h$ is the straight-line distance between two data points.

[0125] Furthermore, the attribute values in this embodiment include effective energy signals, gamma rays, and coal-rock dielectric constants. That is, through the acquisition formula of the semivariogram value, the semivariogram value corresponding to the effective energy signal, the semivariogram value corresponding to the gamma ray, and the semivariogram value corresponding to the coal-rock dielectric constant can be obtained.

[0126] In the process of specific implementation, the first interval grouping contains three data points, and the energy differences between these three data points are: 0.7 J, 0.2 J, and 0.5 J. Therefore, the semivariogram value corresponding to this interval grouping is 0.13 J. 2 。

[0127] Step SA3.2: Obtain the interpolation weights. That is, according to the position of the target voxel center point and the attribute values of all data points within the preset range, establish a Kriging equation system to obtain the weight coefficient corresponding to each data point. At the same time, according to the size of the attribute value corresponding to each data point, determine the attribute interpolation result of the target voxel center point. Specifically as follows:

[0128] Step SA3.2.1: Construct the Kriging equation system. That is, according to the position of the target voxel center point and the positions of all data points within the preset range at the position of the target voxel center point, establish a Kriging equation system through the semivariogram function model to determine the weight coefficient corresponding to each data point.

[0129] In this embodiment, the acquisition formula of the Kriging equation system is specifically:

[0130]

[0131] Where: $\omega$ d is the $d$-th weight coefficient, $\gamma(h$ dp ) is the semivariogram function value between the $d$-th neighboring data point and the $p$-th neighboring data point, $\mu$ is the Lagrange multiplier, $\gamma(h$ d0 ) is the semivariogram function value between the $d$-th neighboring data point and the target voxel center point, $d$ is the index of the weight coefficient and the neighboring data point, $h$ dp is the straight-line distance between the $d$-th neighboring data point and the $p$-th neighboring data point, $h$ d0 is the straight-line distance between the $d$-th neighboring data point and the target voxel center point, and $n$ is the number of neighboring data points participating in the interpolation.

[0132] During the specific implementation process, the position coordinates of the center point of the target voxel are (5.0, 3.0, 2.0), and the preset range is set to 3m. That is, within a 3m range centered on the target voxel center point, 5 data points are obtained. At the same time, the straight-line distances and microseismic energies between each data point and the target voxel center point are respectively: (0.0, 8.2), (0.1, 8.5), (0.55, 7.5), (0.63, 7.8), and (1.2, 8.0).

[0133] Furthermore, according to the obtained straight-line distances and microseismic energies, a Kriging equation system is established, and five weight coefficients can be obtained, and their magnitudes are respectively: 0.3, 0.2, 0.25, 0.15, and 0.1.

[0134] Step SA3.2.2: Obtain the attribute interpolation result. That is, according to the weight coefficients obtained in step SA3.2.1 and the attribute values (effective energy signal, gamma ray, and coal-rock dielectric constant) corresponding to each weight, the attribute interpolation result of the target voxel center point is determined, specifically:

[0135]

[0136] Where: R is the attribute interpolation result of the target voxel center point, n is the number of neighboring data points participating in the interpolation, d is the index of the weight coefficient and the neighboring data point, ω d is the d-th weight coefficient, and Z d is the attribute value of the d-th neighboring data point.

[0137] During the specific implementation process, according to the determined weight coefficients and the microseismic energies corresponding to the above five data points, the corresponding microseismic energy interpolation result is: 7.94J.

[0138] Step SA3.3: Multimodal feature fusion. That is, according to the attribute interpolation results obtained in step SA3.2.2, after normalization processing, the normalized values of the attribute values are obtained. That is to say, the attribute interpolation result of the effective energy signal is normalized to obtain the normalized value of the effective energy signal. Similarly, the attribute interpolation result of the gamma ray is normalized to obtain the normalized value of the gamma ray. Similarly, the attribute interpolation result of the coal-rock dielectric constant is normalized to obtain the normalized value of the coal-rock dielectric constant.

[0139] Furthermore, according to the normalized value of the effective energy signal, the normalized value of the gamma ray, and the normalized value of the coal-rock dielectric constant, the fusion feature value is obtained, specifically:

[0140] L = α1·E nor +α2·ζ nor +α3·εnor

[0141] Among them: L is the fusion eigenvalue, α1 is the weight coefficient of the effective energy signal, and E nor is the normalized value of the effective energy signal, α2 is the weight coefficient of gamma rays, and ζ nor is the normalized value of gamma rays, α3 is the weight coefficient of the coal-rock dielectric constant, and ε nor is the normalized value of the coal-rock dielectric constant.

[0142] In the process of specific implementation, the normalized value of the effective energy signal is 0.082, the normalized value of gamma rays is 1.2, the normalized value of the coal-rock dielectric constant is 0.045, and at the same time, the weight coefficients of the effective energy signal, gamma rays, and the coal-rock dielectric constant are 0.4, 0.3, and 0.3 respectively, then the fusion eigenvalue is 0.4363.

[0143] In this embodiment, the physical modeling module establishes a wave equation regularization layer according to the fusion eigenvalue obtained by the multi-modal perception module, and adjusts the trend of the established joint learning model through the wave equation regularization layer. At the same time, according to the adjusted trend, the final joint learning model is obtained. Specifically as follows:

[0144] Step SB1: Obtain the residual. That is, through the fusion eigenvalue obtained in step SA3.3, obtain the predicted wave velocity, and through the predicted wave velocity and the predicted pressure obtained by the neural network model, obtain the residual between the predicted pressure and the theoretical pressure obtained by the wave equation. Specifically:

[0145]

[0146] Among them: R w is the wave equation residual, p pred is the predicted pressure field corresponding to the neural network model, T is the discretized time, v pred is the predicted wave velocity, e, f, g are three-dimensional spatial grid indices, h is the index of discretized time, is the Laplace operator.

[0147] In this embodiment, the formula for obtaining the predicted wave velocity is specifically:

[0148]

[0149] Among them: v pred is the predicted wave velocity, μ is the Lagrange multiplier, and L is the fusion eigenvalue.

[0150] Step SB2: Obtain the model trend. That is, according to the residual of the wave equation obtained in Step SB1, obtain the trend of the joint learning model, and continuously adjust the trend of the joint learning model based on the obtained model trend. In this embodiment, the formula for obtaining the model trend is specifically:

[0151]

[0152] Where: R w is the residual of the wave equation, p pred is the predicted pressure field corresponding to the neural network model, T is the discretized time, v pred is the predicted wave velocity, is the Laplace operator, and θ is the model trend.

[0153] In this embodiment, the safety decision-making module obtains the final joint learning model through Step SB2, and obtains the risk probability corresponding to each real-time parameter data through the real-time obtained parameter data, and triggers a warning signal according to the obtained risk probability and the warning risk threshold. Specifically as follows:

[0154] Step SC1: Obtain the risk threshold. That is, by the microseismic energy, gas concentration, and stress field concentration obtained in real time, and taking them as the input of the final joint learning model, the corresponding risk probability can be obtained, specifically:

[0155]

[0156] Where: P risk is the risk probability, β1 is the weight coefficient of the signal energy, β2 is the weight coefficient of the gas concentration, β3 is the weight coefficient of the stress gradient, b is the bias term, E is the signal energy of the microseismic event, C is the gas concentration, is the stress gradient.

[0157] Step SC2: Trigger a warning signal. That is, compare the risk probability obtained in Step SC1 with the preset risk threshold, and trigger a warning signal according to the comparison result, specifically:

[0158] When the risk probability is greater than the first-level risk threshold but less than the second-level risk threshold, trigger a first-level warning signal. When the risk probability is not less than the second-level risk threshold, trigger a second-level warning signal. Otherwise, do not trigger a warning signal.

[0159] Although the embodiments of the present invention have been shown and described, for those of ordinary skill in the art, it can be understood that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended embodiments and their equivalents.

Claims

1. A coal mine geological modeling system based on geophysical exploration data processing and machine learning, characterized in that It includes: A multi-modal perception module, which obtains the difference magnitude corresponding to each piece of the multi-source heterogeneous data by collecting the multi-source heterogeneous data, and fuses the multi-source heterogeneous data through the difference magnitude to determine a fused feature value; SA1: Obtain multi-source heterogeneous data streams: Obtain multi-source heterogeneous data within the advancing direction of the excavation face through a microseismic sensor, a drill pipe device, and an electromagnetic wave CT pair; SA2: Perform dynamic time warping: Obtain the coordinates and time corresponding to the effective energy through a microseismic sensor, obtain the interface position and measurement time of the well logging through the drill pipe device, and align the coordinates corresponding to the effective energy and the interface position, and the time corresponding to the effective energy and the measurement time; Step SA3: Obtain a three-dimensional fused feature volume: Obtain corresponding interpolation data through the effective energy obtained by the microseismic sensor, the gamma ray intensity while drilling, and the CT dielectric constant of the electromagnetic wave, and obtain a fused normalized value according to the interpolation data; A physical modeling module, which obtains the pressure residual corresponding to the regularization layer of the wave equation through the fused feature value, and adjusts the trend of the established joint learning model according to the pressure residual to obtain a final joint learning model; A safety decision-making module, which obtains the risk probability corresponding to each piece of the real-time data through the final joint learning model and the real-time data, and triggers an early warning signal according to the risk probability and the early warning risk threshold.

2. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 1, characterized in that, Obtaining multi-source heterogeneous data within the advancing direction of the excavation face includes: SA1.1: Deploy microseismic monitoring nodes: Obtain an effective energy signal through a microseismic sensor and a preset energy threshold, specifically: Where: E is the signal energy of the microseismic event, t1 is the starting point of the integration time window, t2 is the ending point of the integration time window, a x is the instantaneous acceleration value of the microseismic sensor in the x orthogonal direction, a y is the instantaneous acceleration value of the microseismic sensor in the y orthogonal direction, a z is the instantaneous acceleration value of the microseismic sensor in the z orthogonal direction; SA1.2: Set up a logging-while-drilling unit: Correct the azimuth deviation of the detector through the attitude angle of the drill bit and the coordinate transformation matrix, specifically: Where: x′ is the roadway advancing direction, y′ is the transverse direction perpendicular to the roadway, z′ is the gravitational direction vertically downward, R z (ψ) is the yaw angle rotation matrix about the z-axis, ψ is the yaw angle of the drill bit, R y (θ) is the pitch angle rotation matrix about the y-axis, θ is the pitch angle of the drill bit, x is the drilling direction along the axial direction of the drill pipe, y is the transverse direction perpendicular to the axial direction of the drill pipe, and z is the vertical direction perpendicular to the axial direction of the drill pipe; SA1.3: Configure an electromagnetic wave CT pair: Obtain the true coal-rock dielectric constant through an electromagnetic wave transmitter and an array receiver, specifically: Where: ε r is the true dielectric constant of coal and rock, E m is the measured electromagnetic field intensity, E s is the scattered field intensity caused by the metal support structure, E t is the theoretical field strength without metal interference, ||·|| 2 is the square of the norm.

3. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 1, characterized in that, Aligning the coordinates corresponding to the effective energy and the interface position, and the time corresponding to the effective energy and the measurement time includes: SA2.1: Align data timestamps: Obtain the total cost function of time according to the time series of the effective energy signal and the time series of the well logging, specifically: Where: Tot(i, j) is the total alignment cost, Tim(i, j) is the time difference cost, λ is the weight coefficient, IMU(i, j) is the attitude difference cost, i is the time series index of the microseismic event, and j is the time series index of the logging-while-drilling; SA2.2: Spatial mapping: Obtain three-dimensional coordinates in the set roadway coordinate system through UWB tags and mobile tags, and adjust the position of each USB base station according to the three-dimensional coordinates.

4. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 3, wherein Adjusting the position of each USB base station includes: SA2.2.1: Set up a roadway coordinate system: Use the position of the cutting head of the roadheader as the origin of the roadway coordinate system through a UWB tag, use the advancing direction of the roadheader as the X axis, the vertical upward direction of the roadheader as the Y axis, and the horizontal lateral direction of the roadheader as the Z axis to establish a roadway coordinate system; SA2.2.2: Correct the position of the UWB base station: Set UWB base stations at intervals in the roadway, set mobile tags in the microseismic sensors and the drill - while - drilling equipment. At the same time, in the roadway coordinate system, obtain the three - dimensional coordinates of each mobile tag, obtain the Euclidean distance between the three - dimensional coordinates and the UWB base station, and adjust the position of the UWB base station according to the magnitude relationship between the Euclidean distance and the preset error. Specifically: When the Euclidean distance is not less than the preset error, adjust the position of the UWB base station until the Euclidean distance is less than the preset error. Otherwise, do not adjust the position of the UWB base station.

5. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 1, characterized in that, Obtain the fused normalized value, including: SA3.1: Obtain Kriging interpolation: Through the three - dimensional coordinates of each data point, perform interval grouping, and determine the average semivariogram value corresponding to each interval grouping according to the attribute differences of all data points in each interval grouping. SA3.2: Obtain interpolation weights: According to the position of the center point of the target voxel and the attribute values of all data points within the preset range, obtain the weight coefficient of each data point, and determine the attribute interpolation result of the center point of the target voxel according to the weight coefficient. SA3.3: Multimodal feature fusion: Perform normalization processing on the attribute interpolation result, and obtain the fusion feature value according to the normalized attribute interpolation result. Specifically: L = α1·E nor + α2·ζ nor + α3·ε nor Where: L is the fusion eigenvalue, α1 is the weight coefficient of the effective energy signal, E nor is the normalized value of the effective energy signal, α2 is the weight coefficient of the gamma ray, ζ nor is the normalized value of the gamma ray, α3 is the weight coefficient of the coal-rock dielectric constant, ε nor is the normalized value of the coal-rock dielectric constant.

6. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 5, wherein Determine the average semivariogram value corresponding to each interval grouping, including: SA3.1.1: Perform semivariogram function modeling: According to the three - dimensional coordinates of each data point in the roadway coordinate system, obtain the three - dimensional Euclidean distance between data points. Specifically: Where: h is the straight-line distance between two data points, x m is the X-axis coordinate of the m-th data point in the roadway coordinate system, x n is the X-axis coordinate of the n-th data point in the roadway coordinate system, y m is the Y-axis coordinate of the m-th data point in the roadway coordinate system, y n is the Y-axis coordinate of the n-th data point in the roadway coordinate system, z m is the Z-axis coordinate of the m-th data point in the roadway coordinate system, z n is the Z-axis coordinate of the n-th data point in the roadway coordinate system; SA3.1.2: Perform distance grouping: According to the three - dimensional Euclidean distance between data points and the preset distance, perform interval grouping. SA3.1.3: Obtain the average semivariogram value: According to the attribute differences between data points in each interval grouping, obtain the average semivariogram value corresponding to each interval grouping. Specifically: Among them: γ(h) is the semi-variogram value, N(h) is the number of pairs of data points near the lag distance h, Z(X r ) is the attribute value corresponding to the position X r , Z(X r +h) is the attribute value corresponding to another point position h away from the position X r , X r is the three-dimensional coordinate of the r-th data point in the interval grouping in the roadway coordinate system, r is the index of the data point in the interval grouping, and h is the straight-line distance between two data points.

7. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 5, characterized in that Determine the attribute interpolation result of the center point of the target voxel, including: SA3.2.1: Construct the Kriging equation system: According to the position of the center point of the target voxel and all data points within the preset range, establish a Kriging equation system through the semivariogram function model to determine the weight coefficient corresponding to each data point. The Kriging equation system is specifically: where: ω d is the d-th weight coefficient, γ(h dp ) is the semivariogram value between the d-th neighboring data point and the p-th neighboring data point, μ is the Lagrange multiplier, γ(h d0 ) is the semivariogram value between the d-th neighboring data point and the center point of the target voxel, d is the index of the weight coefficient and the neighboring data point, h dp is the straight-line distance between the d-th neighboring data point and the p-th neighboring data point, h d0 is the straight-line distance between the d-th neighboring data point and the center point of the target voxel, and n is the number of neighboring data points participating in the interpolation; SA3.2.2: Obtain the attribute interpolation result: According to the weight coefficient and attribute value corresponding to the data point, determine the attribute interpolation result of the center point of the target voxel. Specifically: Where: R is the attribute interpolation result of the center point of the target voxel, n is the number of neighboring data points participating in the interpolation, d is the index of the weight coefficient and the neighboring data points, ω d is the d-th weight coefficient, Z d is the attribute value of the d-th neighboring data point.

8. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 1, characterized in that, Obtain the final joint learning model, including: SB1: Obtain the residual: Through the fusion feature value, obtain the predicted wave velocity, and obtain the wave equation residual according to the predicted wave velocity and the predicted pressure obtained by the neural network model. Specifically: Where: R w is the residual of the wave equation, p pred is the predicted pressure field corresponding to the neural network model, T is the discretized time, v pred is the predicted wave velocity, e, f, g are three-dimensional spatial grid indices, h is the index of the discretized time, is the Laplace operator; SB2: Obtain the model trend: According to the wave equation residual, obtain the trend of the joint learning model, and adjust the joint learning model according to the trend. The formula for obtaining the trend is specifically: Where: R w is the residual of the wave equation, p pred is the predicted pressure field corresponding to the neural network model, T is the discretized time, v pred is the predicted wave velocity, is the Laplace operator, and θ is the model trend.

9. The coal mine geological modeling system based on geophysical exploration data processing and machine learning according to claim 8, wherein, The formula for obtaining the predicted wave velocity is specifically: Where: v pred is the predicted wave velocity, μ is the Lagrange multiplier, and L is the fusion eigenvalue.

10. The coal mine geological modeling system based on geophysical prospecting data processing and machine learning according to claim 1, wherein Trigger the warning signal, including: SC1: Obtain the risk threshold: Obtain the risk probability through microseismic energy, gas concentration, stress field concentration, and the final joint learning model, specifically: Where: P risk is the risk probability, β1 is the weight coefficient of the signal energy, β2 is the weight coefficient of the gas concentration, β3 is the weight coefficient of the stress gradient, b is the bias term, E is the signal energy of the microseismic event, C is the gas concentration, is the stress gradient; SC2: Trigger the warning signal: Compare the risk probability with the preset risk threshold, and trigger the warning signal according to the comparison result, specifically: When the risk probability is greater than the first-level risk threshold but less than the second-level risk threshold, trigger the first-level warning signal. When the risk probability is not less than the second-level risk threshold, trigger the second-level warning signal. Otherwise, do not trigger the warning signal.

Images