A mine safety production risk prevention and control method based on multi-source data fusion

Abstract

The present application relates to the technical field of safety management, in particular to a mine safety production risk prevention and control method based on multi-source data fusion, comprising the following steps: constructing a coupling parameter set through microseismic and temperature and humidity sensors, extracting displacement points by Bayesian regression detection to draw a deformation map, determining a displacement area, calculating a curvature change rate to locate a mutation area, extracting a thermal strain rate and an axial strain rate to generate a smooth curve to identify a risk point, and redrawing a mine thermal coupling risk map when the overlap rate is more than 60% to expand the boundary. In the present application, a multi-parameter analysis framework is constructed by coupling microseismic and temperature and humidity data to extract a deformation active zone, a difference in air flow velocity gradient is combined with a crack direction to correct an air flow path, a mutation area is located based on a mean threshold of a curvature change rate, a thermal strain and an axial strain are extracted to screen risk points, and a thermal amplitude mapping table is constructed to represent rock mass instability risk.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of safety management, and particularly relates to a mine safety production risk prevention and control method based on multi-source data fusion. BACKGROUND

[0002] Among them, the mine safety production risk prevention and control method based on multi-source data fusion refers to collecting multiple types of information resources such as underground operation environment data, equipment operation data and operation personnel state data, and using time sequence correlation analysis method, statistical feature extraction method and spatial data fusion strategy to integrate and process the same, and then constructing a risk distribution model and early warning information flow of a mine operation area. The method specifically covers dynamic perception of safety critical data such as gas concentration, ventilation volume and support parameters, and uses data fitting function, variable reduction technology and adjacent information correlation reasoning rules to process the data, so as to establish a regional and time period risk control system in mine production.

[0003] In the multi-source data fusion process of the prior art, static time sequence correlation or single statistical feature extraction is used, and the depth of mining of the dynamic coupling relationship of multiple parameters is lacking, including the real-time interactive influence of gas concentration and rock mass displacement, which is not effectively modeled, so that the risk distribution model cannot reflect the dynamic evolution trend under complex environment. The existing method relies on fixed threshold or experience rule for risk warning, and lacks quantitative analysis ability for airflow path mutation and nonlinear characteristics of thermal strain increment trend, and is easy to miss or misjudge, such as when the mine ventilation system is abnormal, the traditional model is difficult to identify the relevance of local airflow disturbance and support structure failure. In addition, the data processing flow of the prior art is isolated, and the environment parameters, equipment state and personnel behavior data do not form a closed loop feedback, and the hidden danger elimination mechanism relies on manual inspection and periodic evaluation, which is not timely, and it is difficult to respond to sudden rock mass deformation or thermal anomaly in time. Including, the coupling effect of sudden increase of gas concentration and rock mass crack expansion is not monitored in real time, which leads to lag in emergency response and cannot start the prevention and control measures in the accident budding stage. The existing spatial data fusion strategy does not fully refer to the unsteady characteristics of the airflow path, and the model boundary condition setting relies on historical data or theoretical assumption, which deviates greatly from the actual dynamic environment, and weakens the reliability of risk warning. SUMMARY

[0004] The purpose of the present application is to solve the shortcomings in the prior art, and a mine safety production risk prevention and control method based on multi-source data fusion is proposed.

[0005] In order to achieve the above purpose, the technical scheme adopted by the present application is as follows: a mine safety production risk prevention and control method based on multi-source data fusion, comprising the following steps:

[0006] S1: Obtain triaxial rock mass displacement rate parameters through microseismic sensors, and collect environmental parameters within the corresponding time through temperature and humidity sensors to construct a coupling parameter set, call a Bayesian linear regression model to perform coupling parameter set anomaly detection and extract displacement points respectively, and construct a rock mass deformation active zone distribution map;

[0007] S2: Determine a displacement area according to the rock mass deformation active zone distribution map, perform difference statistics on the airflow velocity gradient around the displacement area, combine the rock mass displacement rate to extract the rock mass main fracture propagation direction as a boundary condition correction term of the finite volume method model, and construct a non-steady airflow path curve graph;

[0008] S3: Calculate the curvature change rate of the airflow path curve in the non-steady airflow path curve graph, judge whether the curvature change rate of the continuous three sections exceeds 1.5 times of the average value, and locate the local mutation area of the airflow path according to the judgment result;

[0009] S4: Obtain the temperature gradient thermal strain rate and axial strain rate of the rock mass material in the local mutation area of the airflow path, generate the corresponding time sequence smooth curve by using the moving average method, identify the rock mass risk coordinate point group in which the continuous period increments are all positive values, and construct a local thermal amplitude trend graph;

[0010] S5: Call the risk point group of the local thermal amplitude trend graph, calculate the spatial overlap area ratio with the boundary curve of the rock mass deformation active zone distribution map, and if the overlap area ratio is greater than 60%, use a cubic spline interpolation function to fit the boundary curve extension, and redraw a new area to generate a mine thermal coupling risk map.

[0011] As a further scheme of the present application, the rock mass deformation active zone distribution map includes displacement rate anomaly points, environmental parameters, and coupling coefficients, the non-steady airflow path curve graph specifically refers to a velocity gradient difference value, a rock mass main fracture propagation angle, and a corrected boundary condition parameter, the local mutation area of the airflow path specifically refers to a curvature change rate and an airflow velocity anomaly area, the local thermal amplitude trend graph includes temperature strain smooth values, axial strain increment sequences, and risk point spatial coordinates, and the mine thermal coupling risk map includes a spatial overlap area ratio, a cubic spline interpolation boundary parameter, an extended area coordinate set, and a coupling risk level partition.

[0012] As a further scheme of the present application, the specific steps of obtaining triaxial rock mass displacement rate parameters through microseismic sensors, and collecting environmental parameters within the corresponding time through temperature and humidity sensors to construct a coupling parameter set, calling a Bayesian linear regression model to perform coupling parameter set anomaly detection and extracting displacement points respectively, and constructing a rock mass deformation active zone distribution map include:

[0013] S101: Obtain the triaxial displacement rate data of the microseismic sensor, the airflow velocity, temperature gradient and humidity saturation parameters collected by the temperature and humidity sensor, and process the parameters of different dimensions by using the minimum-maximum normalization method, match the four types of parameters based on the unified timestamp, fill in the missing values by using the linear interpolation method, and generate a coupled parameter set;

[0014] S102: Obtain the displacement rate, airflow velocity, temperature gradient and humidity saturation parameters of the coupled parameter set, construct a linear regression equation by using a Bayesian linear regression model, calculate the absolute value of the residual and the 95% confidence interval, screen the measurement point coordinates whose residuals exceed the upper limit of the confidence interval, and generate an abnormal detection result set;

[0015] The Bayesian linear regression model sets the prior distribution of the regression coefficient and the noise, combines the observation data to deduce the posterior distribution, generates the predicted value and the confidence interval of the displacement rate, and is used for residual analysis and abnormal point screening;

[0016] S103: Call the abnormal detection result set, extract the measurement point coordinates with a displacement rate of greater than or equal to 0.5 mm / h, perform density clustering on the measurement points by using a spatial neighborhood clustering algorithm, merge adjacent points with a spacing less than a neighborhood radius, the neighborhood radius is determined by using a K-nearest neighbor method according to the average distance from the measurement point to the Kth nearest neighbor point, and generate a rock mass deformation active zone distribution map.

[0017] As a further scheme of the application, according to the rock mass deformation active zone distribution map, a displacement region is determined, the airflow velocity gradient around the displacement region is differentially counted, the main crack propagation direction of the rock mass is extracted as a boundary condition correction term of a finite volume method model in combination with the rock mass displacement rate, and the specific steps for constructing a non-steady airflow path curve diagram include:

[0018] S201: Call the rock mass deformation active zone distribution map, extract the measurement point coordinates around the displacement region, calculate the airflow velocity gradient based on the Euclidean distance between the measurement points and the airflow velocity difference between adjacent measurement points, and generate a gradient statistical result;

[0019] S202: Call the gradient statistical result, compare the airflow velocity gradient of each measurement point with the airflow velocity gradient perturbation error critical threshold item by item, screen the measurement point coordinates whose airflow velocity gradient exceeds the airflow velocity gradient perturbation error critical threshold, and generate an over-limit region coordinate set;

[0020] The airflow velocity gradient perturbation error critical threshold is determined by parameter sensitivity analysis in the finite volume method, when the gradient exceeds 0.3, and the perturbation error of the airflow to the displacement direction of the rock mass exceeds 5%, the target is taken as a trigger point for boundary condition correction;

[0021] The airflow velocity gradient difference statistics are determined by calculating the Euclidean distance mean of the first-order derivative of the airflow velocity vector in the neighborhood.

[0022] S203: calling the hyperlimit area coordinate set, using a geometric direction statistics method to calculate a feature vector of the rock mass displacement rate, extracting a feature vector direction corresponding to a maximum feature value, taking a main fissure expansion direction of the rock mass as a normal correction quantity of a boundary condition of the finite volume method, and generating a non-steady airflow path curve graph;

[0023] The finite volume method correction term is a projection component of the displacement rate along the main fissure expansion direction of the rock mass superimposed on the original normal velocity component.

[0024] As a further scheme of the present application, the curvature change rate of the airflow path curve in the non-steady airflow path curve graph is calculated, and it is judged whether the curvature change rates of three continuous segments exceed 1.5 times of the average value, and the specific steps of positioning the local mutation area of the airflow path according to the judgment result include:

[0025] S301: obtaining a discrete coordinate point set of the non-steady airflow path curve graph, using a cubic spline interpolation method to parameterize the curve, the condition of the cubic spline interpolation being that the first derivative and the second derivative of the curve at each interpolation node are continuous, using a parameterized cubic spline equation to calculate the curvature value of the non-steady airflow path curve graph point by point, generating a curvature value sequence, and performing difference operation on the curvatures of adjacent sampling points to generate a curvature change rate sequence;

[0026] S302: calling the curvature change rate sequence, calculating the average value of all data points in the curvature change rate sequence, comparing the values of three adjacent data points in sequence with the average value, and if the absolute values of the three data points all exceed 1.5 times of the average value, marking the interval to generate an over-average value marking sequence;

[0027] The curvature change rate is calculated by using a first-order difference method, that is, the curvatures of adjacent sampling points are differentially processed, and the difference value between the curvatures of the current sampling point and the previous sampling point is used to approximately represent the curvature change rate;

[0028] S303: based on the over-average value marking sequence, extracting the starting and ending coordinate point indexes of the marked interval, merging the intervals with an adjacent interval less than five sampling points, determining the boundary range of the curvature mutation area in the airflow path, and positioning the local mutation area of the airflow path.

[0029] Compared with the prior art, the present application has the following advantages and positive effects:

[0030] In the present application, by dynamically coupling the microseismic sensor with the temperature and humidity sensor data, a multi-parameter correlation analysis framework is constructed, the rock mass displacement rate is directly correlated with the environmental parameter change, the nonlinear relationship between parameters is captured by using the Bayesian regression model, the detection sensitivity of abnormal displacement points is improved, and the spatial distribution of the deformation active zone is identified. By differentiating the air flow velocity gradient around the displacement area, combining with the crack propagation direction to correct the boundary conditions of the air flow path model, the dynamic interaction between rock mass displacement and air flow is taken into account, and the physical relevance of the non-steady air flow path prediction is enhanced. The local mutation region is judged based on the mean threshold of curvature change rate, the mutation intensity of the geometric characteristics of the air flow path is quantified, the error of subjective experience interpretation is avoided, and the risk air flow disturbance range is located. The temperature gradient thermal strain and axial strain time series data of the mutation region are extracted, the moving average method is used to suppress noise interference, the risk point group is screened through continuous period increment positive trend, the mapping relationship between thermal amplitude and spatial coordinates is constructed, and the dynamic visualization representation of rock mass thermal instability risk is carried out. The overall scheme forms a full-process closed loop from deformation monitoring to thermal risk prediction through multi-source parameter coupling, dynamic boundary correction, quantitative threshold judgment and trend correlation analysis, and improves the timeliness and accuracy of rock mass instability warning in complex environment. BRIEF DESCRIPTION OF DRAWINGS

[0031] Figure 1 It is a main step schematic diagram of the present application. DETAILED DESCRIPTION

[0032] In order to make the purpose, technical scheme and advantages of the present application more clear, the present application will be further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the present application, and are not used to limit the present application.

[0033] In the description of the present application, it should be understood that the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer" and the like indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present application and simplifying the description, and therefore cannot be understood as indicating or implying that the devices or elements referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as limiting the present application. In addition, in the description of the present application, the meaning of "a plurality of" is two or more, unless otherwise specifically limited.

[0034] Please refer to Figure 1 The present application provides a technical scheme: a mine safety production risk prevention and control method based on multi-source data fusion, comprising the following steps:

[0035] S1: Obtain triaxial rock mass displacement rate parameters through microseismic sensors and collect environmental parameters within the corresponding time period through temperature and humidity sensors, construct a coupled parameter set, call the Bayesian linear regression model to detect anomalies in the coupled parameter set and extract displacement points respectively, and construct a distribution map of active rock mass deformation zones.

[0036] S2: Determine the displacement region based on the distribution map of the active deformation zone of the rock mass, perform differential statistics on the airflow velocity gradient around the displacement region, extract the main fracture propagation direction of the rock mass in combination with the rock mass displacement rate as the boundary condition correction term of the finite volume method model, and construct the unsteady airflow path curve.

[0037] S3: Calculate the rate of curvature change of the airflow path curve in the unsteady airflow path curve diagram, determine whether the rate of curvature change of three consecutive segments exceeds 1.5 times the average value, and locate the local abrupt change area of ​​the airflow path based on the judgment result;

[0038] S4: Obtain the temperature gradient thermal strain rate and axial strain rate of rock material in the local abrupt change area of ​​airflow path, generate the corresponding time-series smooth curve using the moving average method, identify the rock mass risk coordinate point group in the curve where the continuous period increment is positive, and construct a local thermal amplitude trend map.

[0039] S5: Call the risk point group of the local thermal amplitude trend map and calculate the spatial overlap area ratio with the boundary curve of the rock mass deformation active zone distribution map. If the overlap area ratio is greater than 60%, the boundary curve is fitted by the cubic spline interpolation function and the new area is redrawn to generate the mine thermal coupling risk map.

[0040] The distribution map of active rock mass deformation zones includes displacement rate anomalies, environmental parameters, and coupling coefficients. The unsteady airflow path curves specifically include velocity gradient differences, rock mass main fracture propagation angles, and corrected boundary condition parameters. Local abrupt change areas in airflow paths specifically refer to curvature change rates and airflow velocity anomaly zones. The local thermal amplitude trend map includes temperature-strain smoothing values, axial strain increment sequences, and spatial coordinates of risk points. The mine thermal coupling risk map includes the proportion of spatially overlapping areas, cubic spline interpolation boundary parameters, extended area coordinate sets, and coupling risk level zoning.

[0041] Please see Figure 1 This invention provides a technical solution: a method for mine safety production risk prevention and control based on multi-source data fusion, comprising the following steps:

[0042] S101: Acquire triaxial displacement rate data from the microseismic sensor and airflow velocity, temperature gradient, and humidity saturation parameters collected by the temperature and humidity sensor. Process parameters of different dimensions using the minimum-maximum normalization method, perform spatial coordinate matching on the four types of parameters based on a unified timestamp, fill in missing values ​​using linear interpolation, and generate a set of coupled parameters.

[0043] The three-axis displacement rate data of the microseismic sensor and the airflow velocity, temperature gradient, and humidity saturation parameters collected by the temperature and humidity sensor are obtained. The timestamp, three-axis displacement rate (X, Y, Z direction), airflow velocity (m / s), temperature gradient (℃ / m), and humidity saturation (%) parameters are extracted from the raw data collected by the sensor. The X, Y, and Z displacement rates are taken as the arithmetic mean to obtain the comprehensive displacement rate. The maximum values of the four parameters are 2.0 mm / h, 5.0 m / s, 8.0 ℃ / m, and 90%, respectively, and the minimum values are 0 mm / h, 0.1 m / s, 0.5 ℃ / m, and 30%, respectively. The parameters are normalized by calculating the normalized value as (current value - minimum value) / (maximum value - minimum value). The single measurement point displacement rate is 1.2 mm / h, which is normalized as (1.2-0) / (2.0-0) = 0.6. For the measurement points with missing timestamps, the parameter values of the adjacent time points (e.g., 09:00 and 09:05) are extracted and linear interpolation is used to fill in the missing values. The airflow velocity at 09:03 is (3.0 m / s + 5.0 m / s) / 2 = 4.0 m / s. The matched parameters are sorted by timestamp to generate the coupling parameter set, including timestamp, normalized displacement rate, normalized airflow velocity, normalized temperature gradient, and normalized humidity saturation.

[0044] Table 1: Sensor raw parameter collection table

[0045] Time stamp Displacement rate (mm / h) Airflow velocity (m / s) Temperature gradient (°C / m) Humidity saturation (%) 09：00 0.8 3.0 2.5 50 09：05 1.2 5.0 4.0 70

[0046] As shown in Table 1, the raw data is normalized to obtain a displacement rate of 0.4 (normalized value), an airflow velocity of 0.58, a temperature gradient of 0.5, and a humidity saturation of 0.5.

[0047] S102: Obtain the displacement rate, airflow velocity, temperature gradient, and humidity saturation parameters of the coupling parameter set, construct a linear regression equation by a Bayesian linear regression model, calculate the absolute value of the residual and the 95% confidence interval, and screen the measurement point coordinates with residual exceeding the upper limit of the confidence interval to generate an abnormal detection result set;

[0048] The Bayesian linear regression model sets the prior distribution of the regression coefficient and noise, combines the observation data to infer the posterior distribution, generates the predicted value and confidence interval of the displacement rate, and is used for residual analysis and abnormal point screening;

[0049] The coupling parameter set is called, and the prior distribution of the regression coefficient in the Bayesian linear regression model is set as a Gaussian distribution (normal distribution) with a mean of 0 and a variance of 1 (i.e. The prior distribution of the noise variance is an inverse gamma distribution with a shape parameter α = 1 and a scale parameter β = 0.1 (i.e. 2~Inv-Gamma(1,0.1)), by updating the posterior distribution using the normalized displacement rate (dependent variable) and normalized airflow velocity, temperature gradient, and humidity saturation (independent variables) from the observed data, the posterior mean of the regression coefficients is obtained as β1=0.2 (airflow velocity), β2=0.15 (temperature gradient), and β3=0.05 (humidity saturation), and the posterior mean of the noise variance is 0.05. The formula for generating the displacement rate prediction value is as follows: Including the normalized airflow velocity x1 = 0.6, temperature gradient x2 = 0.5, and humidity saturation x3 = 0.4 at measuring point A, the predicted value is... The actual displacement rate y = 0.6, and the absolute value of the residual is |0.6 - 0.215| = 0.385. Calculate the absolute value sequence of the residuals for all measuring points (including 0.385, 0.18, 0.22, and 0.15), sort them, and take the 95th quantile (the 95th position value) as 0.35. Filter the measuring points with residuals exceeding 0.35 (including measuring point A with a residual of 0.385) to generate an anomaly detection result set.

[0050] S103: Call the anomaly detection result set, extract the coordinates of measuring points with displacement rates ≥ 0.5 mm / h, use the spatial neighborhood clustering algorithm to perform density clustering on the measuring points, merge adjacent points with a spacing smaller than the neighborhood radius, and use the K-nearest neighbor method to determine the neighborhood radius based on the average distance from the measuring point to the Kth nearest neighbor point, and generate a distribution map of the rock mass deformation active zone.

[0051] Extract measurement points with a concentrated displacement rate ≥ 0.5 mm / h from the anomaly detection results (normalized value corresponds to original value ≥ 1.0 mm / h). Calculate the Euclidean distance between multiple measurement points, setting K = 5. Calculate the distance between each measurement point and its 5th nearest neighbor, taking the average of all distances (0.8 m) as the neighborhood radius. Merge measurement points with a distance less than 0.8 m, including measurement points C (coordinates 10, 20) and D (coordinates 10.5, 20.3) with a distance of [missing value]. Merge them into the same cluster and output a distribution map of active rock deformation zones, including the coordinates of the cluster center and the coverage area.

[0052] Please see Figure 1 This invention provides a technical solution: a method for mine safety production risk prevention and control based on multi-source data fusion, comprising the following steps:

[0053] S201: Call the distribution map of the active deformation zone of the rock mass, extract the coordinates of the measuring points around the displacement area, calculate the airflow velocity gradient based on the Euclidean distance between the measuring points and the airflow velocity difference between adjacent measuring points, and generate gradient statistical results;

[0054] Call rock mass deformation active zone distribution map, first extract the specified displacement area, including the region K, its boundary coordinates for (5, 15) to (25, 35), and then from the inside and periphery of the area within 2 meters range filter all sensor measuring points, get the measuring point coordinate list, for example, including measuring point A (10, 20), measuring point B (11, 21), measuring point C (12, 22), measuring point D (13, 23) series of measuring points and corresponding geographic coordinate data, then, traverse each pair of adjacent measuring points in the list, the adjacent judging standard is that the spatial distance is less than or equal to the preset maximum adjacent distance, set to 5 meters, for the measuring point pair that meets the condition, including measuring point A and measuring point B, calculate the straight line spatial distance between them, that is, the Euclidean distance, calculate as follows: call measuring point A coordinate (x A ,y A ) = (10, 20) and measuring point B coordinate (x B ,y B ) = (11, 21), distance At the same time, extract the airflow velocity readings collected by the two measuring points at the same time or in a very short time window (including time difference less than 1 minute), get measuring point A airflow velocity v A = 3.0 m / s and measuring point B airflow velocity v B = 3.5 m / s, calculate the airflow velocity difference Δv AB = v B -v A = 3.5-3.0 = 0.5 m / s, divide the velocity difference by the Euclidean distance between them to get the airflow velocity gradient between the two measuring points, the calculation process is as follows: Repeat the distance calculation, velocity difference calculation and gradient calculation process for all selected adjacent measuring point pairs, including measuring point C (12, 22) and D (13, 23), the distance is also d CD ≈1.41 m, if the velocities are v C = 4.0 m / s and v D = 4.8 m / s, then the velocity difference Δv CD = 0.8 m / s, Collect and arrange all the calculated gradient values and their corresponding measuring point pair information (including measuring point number or coordinates) to form a data set including multiple gradient records as the gradient statistical result.

[0055] Airflow velocity gradient calculation formula: The formula is used to calculate the airflow velocity gradient between any two adjacent measuring points i and j, where v i and v j are the airflow velocity readings of measuring points i and j respectively, (x i ,y i ) and (xj y j ) are their coordinates, the absolute difference (or directed difference, depending on the gradient direction definition) of the calculated velocities, the denominator calculates the Euclidean distance between the two points, indicating the rate of velocity change per unit distance.

[0056] Table 2: Airflow velocity gradient statistics table for point pairs

[0057] Pair of measuring points Distance (m) Velocity difference (m / s) Gradient (s-1) A-B 1.41 0.5 0.35 C-D 1.41 0.8 0.57

[0058] As shown in Table 2, the calculated distances, velocity differences, and gradient values for some point pairs are shown.

[0059] S202: Call the gradient statistical results, compare each point airflow velocity gradient with the airflow velocity gradient perturbation error critical threshold item by item, screen the point coordinates whose airflow velocity gradient exceeds the airflow velocity gradient perturbation error critical threshold, and generate the over-limit region coordinate set;

[0060] Call the gradient statistical results generated in the previous step, which include a series of point pairs and their corresponding airflow velocity gradient values, including records (point pair A-B, gradient 0.35s -1 ), (point pair C-D, gradient 0.57s -1 ), obtain the preset airflow velocity gradient perturbation error critical threshold, which is determined according to the parameter sensitivity analysis results of the target. The setting basis is: through finite volume method simulation to simulate the influence of airflow field under different gradient conditions on the displacement direction of rock mass, when the simulation shows that the gradient value reaches 0.3s -1 , and the disturbance error of the displacement direction prediction of the rock mass caused by the airflow due to the gradient first exceeds the maximum error limit value allowed by the project (including being set to 5%), then 0.3s -1 is set as the critical threshold value for triggering boundary condition correction. The parameter sensitivity analysis process involves running multiple sets of simulation with different gradient inputs, recording the corresponding perturbation error percentage, drawing a gradient-error relationship curve, and finding the gradient value corresponding to the error reaching 5%. The threshold value here is 0.3s -1 . Compare the numerical value of each gradient value in the gradient statistical result set with the critical threshold value 0.3s -1 item by item, and perform a screening judgment operation. For point pairs with a gradient value greater than 0.3s -1 , it is determined to be over-limit, including the gradient of point pair A-B is 0.35s -1 , because 0.35>0.3, so the point pair belongs to the over-limit case, the gradient of point pair C-D is 0.57s -1 , because 0.57>0.3, so the point pair also belongs to the over-limit case, if there is a point pair E-F with a gradient of 0.25s -1Because 0.25≤0.3, it does not belong to the over-limit, and the screening operation retains all the gradient values exceeding 0.3 s -1 The coordinates of the involved measuring points of the measuring point pair, since A-B and C-D are both over-limit, the coordinates of measuring points A (10, 20), B (11, 21), C (12, 22), and D (13, 23) are extracted, and the coordinates of all the screened over-limit measuring points are collected to form a current coordinate data set, which is the over-limit region coordinate set.

[0061] Screening condition: gradient > airflow velocity gradient perturbation error critical threshold (0.3 s-1), the condition is used to determine whether the calculated gradient value exceeds the threshold 0.3 s-1 determined based on parameter sensitivity analysis. -1 The associated coordinates of the measuring point pair that meets this condition will be screened out.

[0062] S203: Call the over-limit region coordinate set, use the geometric direction statistical method to calculate the feature vector of the rock mass displacement rate, extract the feature vector direction corresponding to the maximum eigenvalue, and take the rock mass main crack propagation direction as the normal correction amount of the finite volume method boundary condition to generate a non-steady airflow path curve graph.

[0063] The finite volume method correction term is the original normal velocity component superimposed with the projection component of the displacement rate along the rock mass main crack propagation direction.

[0064] Call the over-limit region coordinate set generated in the previous step, which includes a series of measuring point coordinates judged to be airflow velocity gradient anomalies, such as measuring points A (10, 20), B (11, 21), C (12, 22), and D (13, 23). Extract the rock mass displacement rate data corresponding to the over-limit measuring points, which include the displacement rate components in the X and Y directions (or the X, Y, and Z directions in three-dimensional space). The displacement rate of measuring point C is (v Cx ,v Cy ) = (1.2 mm / h, 0.8 mm / h), and the displacement rate of measuring point D is (v Dx ,v Dy ) = (1.5 mm / h, 0.9 mm / h). Organize the displacement rate component data of all measuring points in the over-limit region to construct an n x d displacement rate matrix, where n is the number of over-limit measuring points and d is the spatial dimension (2 or 3). Then, based on the displacement rate matrix, calculate the covariance matrix C. For the two-dimensional case, the covariance matrix is a 2 x 2 symmetric matrix, and the elements are calculated as follows: the diagonal elements are the variances of the directional components, and the non-diagonal elements are the covariances between different directional components. Using the data of measuring points C and D (sample size is 2), Cov(X,X) = 0.045, Cov(X,Y) = 0.015, and Cov(Y,Y) = 0.005 are calculated, and the covariance matrix is Then, the covariance matrix C is eigenvalue decomposition, the eigenvalue λ and the corresponding eigenvector v are calculated, the characteristic equation det(C-λI)=0 is solved, two eigenvalues λ1=0.055 and λ2≈-0.00045 (close to 0 in the calculation accuracy) are calculated for the above C, the largest eigenvalue λ max =0.055 is selected, the eigenvector corresponding to the largest eigenvalue λ max is extracted, the eigenvector v1=(k·0.894, k·0.447) is obtained by solving the linear equation group (C-λ max I)v=0, the unit vector (0.894, 0.447) is taken, the direction of the eigenvector represents the direction in which the data changes most in the two-dimensional space, that is, the main direction of the rock mass displacement rate change, that is, the inferred main crack propagation direction, the main crack propagation direction (including the direction vector (0.894, 0.447), and the corresponding azimuth is about 26.6°) is used as a correction amount and applied to the boundary condition setting of the finite volume method, the boundary normal parameter is adjusted to be consistent with or associated with the main crack direction detected, and finally, based on the corrected boundary condition, the finite volume method simulation calculation is re-run to generate a non-steady airflow path curve diagram considering the influence of the rock mass structure.

[0065] Eigenvalue decomposition (main direction determination): Cv=λv, the formula is the standard form of the eigenvalue problem, which is used to solve the eigenvalue λ and the eigenvector v of the covariance matrix C, and the goal is to find the eigenvector v1 corresponding to the largest eigenvalue λ max The direction of the vector indicates the direction in which the data variance is the largest, that is, the key direction of the rock mass displacement.

[0066] Please refer to Figure 1 The application provides a technical scheme: a mine safety production risk prevention and control method based on multi-source data fusion, comprising the following steps:

[0067] S301: Obtain a discrete coordinate point set of a non-steady airflow path curve diagram, parameterize the curve by using a cubic spline interpolation method, the condition for cubic spline interpolation is that the first derivative and the second derivative of the curve at each interpolation node are continuous, the curvature value of the non-steady airflow path curve diagram is calculated point by point by using a parameterized cubic spline equation, a curvature value sequence is generated, and the curvature of adjacent sampling points is subjected to difference operation to generate a curvature change rate sequence;

[0068] The non-steady airflow path curve generated in the preceding step is obtained in the form of a series of discrete two-dimensional or three-dimensional coordinate points, which represent the trajectory of the airflow in space over time. For example, the obtained point set is P0(0, 0), P1(1, 2.5), P2(3, 1.5), P3(5, 3.0), P4(6, 1.0), and P5(7, 2.0). These points form the skeleton of the airflow path. In order to obtain a smoother and differentiable curve representation and calculate its geometric properties, a cubic spline interpolation technique is used to process these discrete points to construct a piecewise cubic polynomial curve S(t) = (x(t), y(t)) that passes through all given points. The parameter t can be taken as the cumulative chord length or directly use the point index as the parameter node t0 = 0, t1 = 1,..., t5 = 5. The core condition for constructing the spline curve is that at each internal node P i i = 1, 2, 3, 4, the values, first derivative values (representing the tangent direction), and second derivative values (related to curvature) of the adjacent two cubic polynomials S i-1 (t) and S i (t) must be equal, that is, S i-1 (t i ) = S i (t i ) = P i , S′ i-1 (t i ) = S′ i (t i ), and S″ i-1 (t i ) = S″ i (t i ). Usually, boundary conditions are also added, including the second derivative at the end points of the natural spline hypothesis curve S″0(t0) = 0 and S″ n-1 (t n ) = 0. Solving the linear equation system formed by these conditions can determine the coefficients of each cubic polynomial, obtaining the parameterized cubic spline equations x(t) and y(t). Then, using the obtained parameterized equations, the curvature κ(t) of each point on the curve is calculated. Curvature is a geometric quantity that measures the degree of bending of the curve, and its calculation requires the first and second derivatives of the parameter equation, which calls x′(t), y′(t), x″(t), and y″(t) to substitute into the curvature calculation expression A series of dense sampling points are selected within the definition domain of the parameter t, for example, from t = 0 to t = 5, with a step size of Δt = 0.1, obtaining t k = k × 0.1. For each t k , the corresponding curvature value κ k = κ(t k ) is calculated, forming a curvature value sequence {κ0, κ1,..., κ50} in order to analyze the curvature variation, a first-order difference operation is performed on this sequence, and the curvature difference between adjacent sampling points, i.e. the curvature variation rate Δκ k = κ k - κ k-1 , is calculated, where k starts from 1 to 50. This difference value approximately represents the rate of curvature change when the parameter t changes from t k-1 to t k . All the calculated difference values Δκ k are collected to generate the final curvature variation rate sequence.

[0069] Curvature calculation formula: This formula calculates the curvature κ(t) of the parameterized planar curve (x(t), y(t)) at the parameter t. The parameters x'(t) and y'(t) are the first-order derivatives of the x and y coordinates of the curve with respect to the parameter t, representing the tangent vector components of the curve at this point. The parameters x"(t) and y"(t) are the second-order derivatives, related to the acceleration or bending manner of the curve. The numerator |x'y"-y'x"| is the two-dimensional form of the cross product of the tangent vector and the acceleration vector (if t is time), and the denominator (x' 2 +y' 2 ) 3 / 2 is the cube of the speed (tangent vector modulus). The entire formula is derived based on differential geometry, and the absolute value of the calculation result κ(t) represents the bending degree of the curve at this point, with a larger value indicating a more severe bend. Parameter assignment and acquisition: Assuming that the derivative values at a point t = 1.5 obtained through cubic spline interpolation are as follows: x'(1.5) = 2.0, y'(1.5) = -0.5, x"(1.5) = 0.8, and y"(1.5) = 1.2. These values are calculated by taking the derivative of the interpolated spline function, which is the derivative of the continuous function obtained after interpolating the original discrete points P0 to P5. Formula calculation: Substitute the above derivative values into the curvature formula: The benefit of this formula is that by accurately calculating the curvature of each point of the curve, the local bending degree of the airflow path can be quantified, providing a basis for subsequent identification of path mutation regions. Result interpretation: The calculated curvature value κ(1.5) ≈ 0.3195 represents the bending degree of the airflow path curve at the parameter t = 1.5. The size of this value itself needs to be compared with the curvatures of other points, and a larger curvature value usually corresponds to a sharp turn on the path. Calculate the curvature values of all sampling points to form a sequence {κ k}, such as {…, 0.15, 0.32, 0.45, 0.20, …}, which is the basis for subsequent calculation of the curvature variation rate.

[0070] Table 3: Parameterized airflow path curvature and curvature variation rate table

[0071] Parameter point t k ]] Curvature K k ]] Rate of change of curvature AK k ]]> 1.4 0.15 - 1.5 0.32 0.17 1.6 0.45 0.13 1.7 0.20 -0.25

[0072] Table 3 shows the curvature of some sampling points and the rate of change of curvature calculated by the first-order difference.

[0073] S302: Call the curvature change rate sequence, calculate the mean of all data points in the curvature change rate sequence, compare the values ​​of three consecutive adjacent data points in the sequence with the mean, and if the absolute values ​​of the three data points all exceed 1.5 times the mean, mark the interval and generate the over-mean marked sequence.

[0074] The rate of change of curvature is calculated using the first-order difference method, which involves differentiating the curvature values ​​of adjacent sampling points and approximating the rate of change of curvature by the difference between the curvature of the current sampling point and the curvature of the previous sampling point.

[0075] Call the curvature change rate sequence generated in the previous step, denoted as {Δκ1,Δκ2,...,Δκ}. M Let M = 50 be the length of the sequence. For example, the sequence is {0.1, -0.2, 0.17, 0.13, -0.25, 0.6, 0.7, 0.8, -0.1, 0.2, ...} (containing 50 values). First, calculate the arithmetic mean of the absolute values ​​of all data points in this sequence, also known as the mean absolute deviation (if the center point is 0) or simply the mean of absolute values. The calculation process is to add the absolute values ​​of each element in the sequence and then divide by the length of the sequence, M. Assuming the calculated AbsMean = 0.30, a judgment threshold is set, which is equal to 1.5 times the calculated mean. The threshold is calculated as Threshold = 1.5 × AbsMean = 1.5 × 0.30 = 0.45. This threshold is set to identify data points whose rate of change deviates significantly from the average level. The factor of 1.5 is an empirical coefficient used to balance the sensitivity and false alarm rate of detection. Then, the curvature change rate sequence is traversed, and a sliding window containing three consecutive adjacent data points is used for inspection. For the k-th point in the sequence (k ranges from 2 to M-1), it and its immediate neighbors are checked, i.e., Δκ. k-1 ,Δκ k ,Δκ k+1 For these three points, take the absolute value of each value, and then compare each absolute value with the calculated threshold of 0.45. The following logic is applied: the condition |Δκ| is satisfied if and only if the absolute values ​​of all three data points are greater than the threshold of 0.45. k-1 |>0.45 and|Δκ k |>0.45 and|Δκ k+1|0.45, the index k corresponding to the center point of the current window and its front and back points k-1, k+1 are marked as belonging to the potential curvature mutation interval, for the sequence {-0.25, 0.6, 0.7, 0.8, -0.1}, check the window (-0.25, 0.6, 0.7) centered at k = 6 (value 0.6): |-0.25| = 0.25≯0.45, does not meet the condition, check the window (0.6, 0.7, 0.8) centered at k = 7 (value 0.7): |0.6| = 0.6>0.45, |0.7| = 0.7>0.45, |0.8| = 0.8>0.45, all three conditions are met, so mark indexes 6, 7, 8, check the window (0.7, 0.8, -0.1) centered at k = 8 (value 0.8): |-0.1| = 0.1≯0.45, does not meet the condition, after traversing the entire sequence, a binary (0 or 1) marking sequence is generated according to the marking result, the length is the same as the original sequence, where the marked index position is 1 and the unmarked is 0, for the above-mentioned segment, the positions corresponding to indexes 6, 7, 8 in the marking sequence are 1, and the rest are 0, this sequence is the super-mean marking sequence.

[0076] Super-mean threshold calculation: This formula calculates the threshold for screening significant curvature change rate, M is the total length of the curvature change rate sequence {Δκ k}, Calculate the sum of the absolute values of all elements in the sequence, Get the average value of the absolute value AbsMean, and finally multiply this average value by the factor 1.5 to get the final threshold Threshold, the selection of the factor 1.5 is based on experience or analysis of noise level and signal strength in a specific application scenario, aiming to identify changes higher than 1.5 times the average fluctuation level. Parameter assignment and acquisition: assume that the curvature change rate sequence {Δκ k} contains M = 50 data points, the average value of its absolute value is AbsMean = 0.30. Formula operation:

[0077] Threshold = 1.5 x 0.30 = 0.45, the advantage of the formula is that: by calculating the mean of the absolute value of the data sequence itself and multiplying a coefficient to set the threshold, so that the threshold can be adaptive to the overall volatility level of the data, rather than a fixed hard threshold, which improves the adaptability to different airflow path characteristics. Result interpretation: the calculated threshold Threshold = 0.45 will be used for subsequent comparison and judgment, using this threshold to screen out the area where the curvature change rate sequence of the three consecutive points is significantly deviated from the average change level. The result shows that only when the absolute value of the curvature change rate of the three consecutive sampling points is greater than 0.45, it is considered that there is a significant curvature mutation at this place. This threshold is the key judgment basis for generating the super-mean mark sequence.

[0078] S303: Based on the super-mean mark sequence, extract the start and end coordinate point index of the marked interval, merge the intervals with adjacent intervals less than five sampling points, determine the boundary range of the curvature mutation area in the airflow path, and locate the local mutation area of the airflow path.

[0079] Based on the super-mean marking sequence generated in the previous step, the sequence is a binary sequence, where the position with a value of 1 indicates that the corresponding point in the original curvature rate sequence satisfies the condition that the absolute values of three consecutive points are all more than 1.5 times the mean value. For example, the obtained marking sequence is {0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0} (the index starts from 1), first, identify and extract all continuous 1 subsequences in the sequence, these subsequences represent the initially identified curvature change regions, for the example sequence, three intervals can be identified: the first interval is from index 4 to 6, the second interval is from index 9 to 10, and the third interval is from index 14 to 17, record the start index and end index of each interval, get the interval list [4, 6], [9, 10], [14, 17], next, check the time or space interval (measured in the number of sampling points) between these identified intervals, calculate the interval length between adjacent two intervals, that is, the difference between their end index and the start index of the next interval minus 1, for example, the interval [4, 6] and the interval [9, 10] have an interval of indexes 7 and 8, a total of (9-6)-1=2 sampling points, the interval [9, 10] and the interval [14, 17] have an interval of indexes 11, 12, 13, a total of (14-10)-1=3 sampling points, set an interval merging threshold, which specifies that if the interval between two adjacent marking intervals is less than 5 sampling points, it is considered that the two intervals actually belong to the same large mutation region and need to be merged, compare each interval length calculated with the threshold 5, for the interval [4, 6] and [9, 10] with an interval of 2, because 2<5, it meets the merging condition, merge them into a new interval, take the start index 4 of the first interval as the start index and take the end index 10 of the second interval as the end index, get the merged interval [4, 10], for the merged interval [4, 10] and the next interval [14, 17], their interval is indexes 11, 12, 13, a total of (14-10)-1=3 sampling points, because 3<5, it also meets the merging condition, take the start index 4 of the previous interval as the start index and take the end index 17 of the next interval as the end index, get the final merged interval [4, 17], repeat this merging process until there is no interval between adjacent intervals less than 5 sampling points, the final merged interval list represents the boundary range of the region where the curvature of the airflow path changes significantly and continuously, for example, the final determined mutation region boundary range is index from 4 to 17, according to these index ranges, locate the corresponding space paragraphs on the original airflow path coordinate point set or parameterized curve, these corresponding space paragraphs are the local mutation regions of the airflow path to be found.

[0080] Interval merging condition: (Index start,i+1-Index end,i )-1<5, this formula is used to determine whether two adjacent, preliminary identified curvature mutation intervals (the i-th interval ends at Index end,i , the i+1-th interval starts at Index start,i+1 ) need to be merged, calculate the number of interval sampling points between them (i.e. the number of intermediate points not containing the interval endpoints), compare this interval number with the preset threshold 5, if the interval is less than 5 (i.e. the interval point number is 0, 1, 2, 3 or 4), then perform the merging operation, the selection of threshold 5 is based on the consideration of the airflow path characteristics and the sampling density, aiming to connect those spatially close enough, possibly belonging to the same physical phenomenon (such as bypassing obstacles or entering narrow channels) caused continuous curvature changes.

[0081] Please refer to Figure 1 , the present application provides a technical solution: a mine safety production risk prevention and control method based on multi-source data fusion, comprising the following steps:

[0082] S401: call airflow path local mutation region coordinates, extract temperature gradient thermal strain rate and axial strain rate original sampling point data of rock mass material based on coordinate range, interpolate and align time series data of two kinds of strain rates, map strain rate values to 0-1 interval using linear normalization method, generate time-aligned strain rate data set;

[0083] Call the spatial coordinate range of the airflow path local mutation region determined in the previous step, including the region R identified, whose spatial boundary is determined by the index range [4, 17], or defined by the specific three-dimensional coordinate box

[0084] [x min ,x max ],[y min ,y max ],[z min ,z max ] defined, for example [8, 15], [18, 25], [5, 10] meters, according to this spatial range, filter and extract the original sampling point data recorded by the sensors installed inside or adjacent to the boundary of the region from the rock mass monitoring database, especially the two kinds of strain rate data related to the mechanical behavior of rock mass material, namely temperature gradient thermal strain rate ∈ th (10 -6 / h) and axial strain rate ∈ ax (10 -6 / h), obtain the time series records of the two groups of data, for example, the temperature gradient thermal strain rate sequence E th =

[0085] {(t th,1 ,∈ th,1 ),(tth,2 ,∈ th,2 ),…}={(0h,1.2),(2h,1.5),(4h,1.4),…}, axial strain rate sequence E ax =(t ax,1 ,∈ ax,1 ),(t ax,2 ,∈ ax,2 ),…}={(0h,3.5),(1h,3.8),(3h,4.0),(5h,4.2),…}, since the sampling time points t th,i and t ax,j of the two groups of data can be different, direct comparison or joint analysis is difficult, and it is necessary to align their time bases, adopt time series interpolation processing, select linear interpolation method, and determine a common, higher frequency time grid, for example, with 0.5 hours interval, T sync =

[0086] {0h,0.5h,1.0h,1.5h,…}, for the temperature gradient thermal strain rate sequence E th , the value at t=1.0h needs to be interpolated, which is between t th,1 =0h and t th,2 =2h, and the linear interpolation calculation is The value of the axial strain rate sequence E ax at t=2.0h also needs to be interpolated, which is between t ax,2 =1h and t ax,3 =3h, and the calculation is For all time points on T sync that need to be interpolated, two time-aligned sequences E′ th ={(t k ,∈′ th,k )} and E′ ax ={(t k ,∈′ ax,k )} are obtained, where all t k ∈T sync , next, in order to eliminate the influence of strain rate dimension and numerical range difference, the aligned two sequences are respectively subjected to linear normalization processing, first, the minimum value ∈′ min and the maximum value ∈′ max of each sequence are found, for example, for the aligned E′ th sequence, assuming that the minimum value is ∈′ th,min =0.8 and the maximum value is ∈′ th,max =2.5, for the aligned E′ ax sequence, the minimum value is ∈′ax,min =3.0, the maximum value is ∈′ ax,max =5.5, then, for each data point ∈′ in the sequence k Apply normalization transformation and calculate the normalized value. For example, for ∈′ th (1.0) = 1.35, normalized value is For ∈′ ax (2.0) = 3.9, normalized value is This normalization calculation is performed on all aligned data points, ultimately generating a time-aligned strain rate dataset containing two normalized time series.

[0087] Linear normalization formula: This formula is used to transform the original (or interpolated) time series data points ∈′ k Mapped to the interval [0, 1], ∈′ min It is the minimum value of the time series within the scope of the study, ∈′ max It corresponds to the maximum value minus the minimum value (∈′). k -∈′ min ) Shift the data to start from 0, then divide by the range of the data (∈′). max -∈′ min ) Scale the data to between 0 and 1, where ∈′ k Corresponding ∈′ th,k or ∈′ ax,k The corresponding ∈′ min and ∈′ max These also correspond to the minimum and maximum values ​​of their respective sequences. Parameter assignment and retrieval: Assume the aligned axial strain rate sequence E′ ax At a certain time point t k The value is ∈′ ax,k =4.5×10 -6 / h, the minimum value of this sequence within the time period under consideration is ∈′ ax,min =3.0×10 -6 / h, the maximum value is ∈′ ax,max =5.5×10 -6 / h. These values ​​are obtained by iterating through the aligned time series E′. ax Obtained. Formula calculation: ∈″ ax,k =0.6. The advantage of this formula is that by mapping strain rate data of physical dimensions or numerical ranges to a unified interval of [0,1], scale differences are eliminated, making subsequent comparisons or joint processing of two sequences more fair and convenient, and preventing the sequence with larger numerical values ​​from dominating the analysis results. Interpretation of results: The calculated normalized value ∈″ ax,k=0.6 indicates the original axial strain rate value is 4.5 × 10⁶. -6 / h occupies the 60% relative position within its entire range of variation (from 3.0 to 5.5). After normalizing the data for all time points, the resulting time-aligned strain rate dataset {(t k ,∈″ th,k ,∈″ ax,k This can be used for subsequent smoothing and trend analysis.

[0088] Table 4: Normalized Strain Rate Data

[0089]

[0090]

[0091] Table 4 shows the original (interpolated) strain rate data at some time points and the normalized values ​​obtained by linear normalization.

[0092] S402: Based on the time-aligned strain rate dataset, establish sliding windows for the temperature gradient thermal strain rate and axial strain rate respectively. Define the window to cover the first two and last two data points. Calculate the arithmetic mean of the five data points in the window to replace the original center point value and generate a smoothed strain rate curve set.

[0093] The window size definition of the moving average method is based on the signal main period (FFT acquisition) as the benchmark, combined with the noise standard deviation evaluation to set the lower limit, and the trend continuity is ensured through autocorrelation verification, and dynamically adjusted to balance noise suppression and trend capture.

[0094] Based on the time-aligned and normalized strain rate dataset generated in the previous step {(t k ,∈″ th,k ,∈″ ax,k )}, respectively for the temperature gradient thermal strain rate sequence {∈″ th,k} and axial strain rate sequence {∈″ ax,k The moving average method is applied by creating a sliding window for each sequence. The window size is set to 5 data points, with the current data point k at the center and including the two preceding data points k-2 and k-1 and the two following data points k+1 and k+2. The arithmetic mean of these five data points is calculated, and this mean is used as the new value of the smoothed sequence at point k. The calculation process is as follows: For example, for a normalized axial strain rate sequence, at time point t3 (index k = 3), its original value is ∈″. ax,3 =0.320, assuming the data points before and after it are respectively ∈″ ax,1 =0.200,∈″ ax,2 =0.260,∈″ ax,4=

[0095] 0.360,∈″ ax,5 = 0.360 (from Table 5 and subsequent assumption), then the smoothed value For the start and end portions of the sequence (e.g. k = 1, 2 and k = M-1, M), since there are not enough 5 points available, the window needs to be adjusted or the boundary needs to be handled, for example, for k = 1, the average of the first 3 points (k = 1, 2, 3) can be used, or for k = 2, the average of the first 4 points (k = 1, 2, 3, 4) can be used, or the data points can be supplemented by methods such as mirror filling and then the 5-point window is applied. The size (width of 5) of the sliding window is not randomly set, and its determination includes: first, the original strain rate signal (before normalization) is analyzed by fast Fourier transform (FFT) to identify the main periodic component of the signal, for example, the main period is about 10 sampling points, and a reference window size is set, for example, half of the main period, i.e. 5 points; then the noise level of the signal is evaluated, the standard deviation σ of the original signal (or its difference) is calculated, for example, σ = 0.05 (normalized unit), and the lower limit of the window size is set, for example, the window size is at least able to cover the fluctuation range of ±1σ, and when the window size is 5, the smoothing effect is equivalent to low-pass filtering, which can suppress high-frequency noise with a period less than 5 sampling points; then the autocorrelation functions of the original sequence and the smoothed sequence are calculated, and their decay rates are compared, if the autocorrelation function after smoothing decays too fast, it means that the window is too large and too much trend information is lost, the window needs to be reduced; the window size can also be dynamically adjusted according to the real-time data characteristics, the window size is appropriately increased when the noise is large, and the window size is reduced when fast changes need to be captured, and here it is fixed at 5. This moving average calculation is applied to all data points of the two normalized sequences to generate two smoothed strain rate curves and

[0096] The moving average calculation formula (window size is 5):

[0097] This formula calculates the 5-point center moving average value of the kth point in the sequence ∈″ k+j is the value of the original (normalized) sequence at index k+j, and the summation symbol indicates that the values of the 5 points from k-2 to k+2 are added up, and then divided by the window size 5 to get the arithmetic mean, this average will replace the original value ∈″ kAs the value of the smoothed sequence at the kth point, this process is performed for each point in the sequence (except for the boundary points), effectively filtering out random noise and making the curve smoother. Parameter assignment and acquisition: using the normalized axial strain rate sequence {∈″ ax,k} obtained in S401, the 5 points around index k = 3 are:

[0098] ∈″ ax,1 = 0.200, ∈″ ax,2 = 0.260, ∈″ ax,3 = 0.320, ∈″ ax,4 = 0.360, ∈″ ax,5 = 0.360. These values come from the calculation results of the previous step. Formula operation: The advantage of the formula is that by calculating the average value of the local data points to replace the center point value, it can effectively suppress the high-frequency noise and random fluctuations in the time series, making the potential trend and periodic changes more clear, providing more reliable basic data for subsequent analysis of trend changes (such as calculating the difference). Result interpretation: the smoothed value is the result of the 5-point moving average of the axial strain rate at index 3, which is slightly lower than the original value of 0.320, reflecting the average effect of the surrounding lower values. By replacing all points in the sequence (after processing the boundaries) with the corresponding moving average values, the smoothed axial strain rate curve and the thermal strain rate curve are obtained, which constitute the smoothed strain rate curve set.

[0099] S403: Based on each time node of the smoothed strain rate curve set, calculate the difference between the temperature gradient thermal strain rate and the adjacent node of the axial strain rate, extract the coordinate points with a difference of greater than zero for three consecutive periods, and project the coordinate points that meet the difference of greater than zero for three consecutive periods to a two-dimensional grid with the horizontal axis as the time axis and the vertical axis as the strain rate amplitude, and construct a local thermal amplitude trend chart.

[0100] Based on the smoothed strain rate curve set generated in the previous step, which includes the smoothed temperature gradient thermal strain rate sequence and the axial strain rate sequence For each time node t k (k from 2 to M, where M is the sequence length) represented on the two smoothed curves, calculate the numerical difference between it and the previous adjacent node t k-1 , i.e. calculate the first-order backward difference, for the temperature gradient thermal strain rate, calculate For the axial strain rate, calculate These two differences respectively represent the change in the temperature gradient thermal strain rate and the axial strain rate in the time interval [t k-1 , t kthe change amount or approximate change rate of the two strain rates, next, a condition is set for screening specific change trends: find the time points at which both strain rates show a continuous growth trend, specifically find the index k in the sequence such that the continuous three differences (corresponding to three consecutive time periods or sampling intervals) from this point are all greater than zero, that is, the judgment condition is: and and and simultaneously satisfy and and For all possible starting indexes k (from 2 to M-2), this condition is judged, and all starting index k values that satisfy this condition are extracted, and the coordinate points (i.e. time points t k ) corresponding to these indexes are identified. Assuming that the index set that satisfies the condition obtained through screening is K inc ={k1,k2,…}, for example K inc ={5,18,35}, these identified coordinate points are projected onto a two-dimensional grid plane according to their spatial positions (although this is mainly dealing with time series, the coordinate points originally mean indexes or time points in the time series), the horizontal axis of the grid is set as the time axis, representing time t k , and the vertical axis is set as the strain rate amplitude, which can be selected as one of the strain rates (such as the axial strain rate ) or some combination (such as their sum or average) of the two strain rates as the vertical coordinate value. For each index k∈K inc that satisfies the condition, a point is plotted on the two-dimensional grid with coordinates (assuming the axial strain rate amplitude is selected), all these projected points are connected or displayed in a scatter plot, and finally a graphical representation is constructed, i.e. a local thermal amplitude trend graph, which visually displays the time points at which both strain rates in the identified airflow path mutation region show a continuous growth trend and their corresponding strain rate levels.

[0101] Continuous growth judgment condition: This logical expression is used to determine whether a trend of continuous growth of both strain rates for at least three sampling periods has started at time index k, and are the first-order differences of the two smoothed strain rates, ∧ represents logical AND (AND), and the entire condition requires that the difference values of the two strain rates in the three time intervals from index k to k+2 must all be positive, and the difference values greater than zero indicate that the strain rate is increasing in that time interval.

[0102] Please refer to Figure 1The application provides a technical scheme: a mine safety production risk prevention and control method based on multi-source data fusion, comprising the following steps:

[0103] S501: Obtain the risk point group of the local thermal amplitude trend map and the boundary curve of the rock mass deformation active zone distribution map, convert the vector boundary into raster data, count the number of pixels of the intersection of the raster where the risk point group is located and the boundary raster, calculate the proportion of the intersection pixel number in the total risk point group pixel number, and generate the overlap area proportion;

[0104] Call the risk point group data of the local thermal amplitude trend map and the vector boundary curve of the rock mass deformation active zone distribution map, the risk point group data coordinates are discrete point set {(x i ,y i )}, including points (120, 85), (125, 90), (130, 88), …, and the vector boundary curve of the rock mass deformation active zone is composed of node coordinates of multiple lines, for example, the node sequence {(100, 80), (110, 85), (120, 90), (135, 95), (150, 85)}, convert the vector boundary curve into raster data, set the raster resolution to 1 meter, generate a raster matrix covering the entire analysis area, and the size of each raster unit is 1 meter*1 meter; the raster attribute value is 0 (non-boundary) or 1 (boundary); the specific operation is to traverse all raster units, judge whether the center point coordinates are located in the polygon coverage area of the vector boundary curve, for example, the raster unit (x g ,y g )=(121, 86), the center point coordinates are (121.5, 86.5), if the point falls within the vector boundary polygon, the raster attribute value is set to 1, otherwise, it is 0, and the risk point group coordinates are converted into raster data, each risk point corresponds to the raster unit where it is located, and the attribute value is accumulated; the point (120, 85) is located in the raster unit (120, 85), and the attribute value is added by 1; the number of pixels of the intersection of the raster where the risk point group is located and the boundary raster is counted, that is, the total number of raster units that simultaneously satisfy the risk point attribute value greater than 0 and the boundary attribute value of 1, for example, the total number of risk point group rasters is 200, and the number of units overlapping with the boundary raster is 130; the overlap area proportion is calculated

[0105] Raster intersection calculation: traverse each raster unit (x g ,y g ), judge whether the following conditions are met simultaneously: 1. The risk point group raster attribute value V risk (x g ,y g )>0 2. The boundary raster attribute value V boundary (x g ,y g )=1 If satisfied, the counter N overlapAdd 1, final calculation where N risk is the total number of risk point group grids.

[0106] Parameter assignment and acquisition: the total number of risk point group grids is 200, the total number of boundary grid units is 180, and the number of overlapping units is calculated by traversal to be 130.

[0107] Threshold setting basis: the critical threshold of 60% is set by analyzing mine disaster data, including statistics of 20 rock mass instability events, it is found that when the overlap area ratio R overlap of risk points and deformation boundaries is greater than or equal to 60%, the corresponding rock failure probability P failure reaches 92% (18 / 20 events), while when R overlap <60%, the failure probability is only 15% (3 / 20 events), the threshold of 60% is taken as the critical point of high probability failure, and its significance (p<0.01) is verified by a logistic regression model, and the formula is: where β0=-5.2, β1=8.3, when R overlap =0.6: In actual engineering, it is adjusted to 60% to ensure safety redundancy.

[0108] P failure is the rock failure probability, R overlap is the overlap area ratio of risk points and boundaries, which represents the statistical results of disaster data, β0 is the logarithmic ratio of the baseline failure probability, β1 is the influence coefficient of the overlap area on the failure probability, and the data fitting result reflects the increase of the failure probability logarithmic ratio per 1% increase of the overlap area.

[0109] Formula calculation: Result interpretation: R overlap =0.65 exceeds the critical threshold of 0.6, triggering the boundary expansion operation.

[0110] Table 5: Risk point group and deformation boundary grid overlap statistics

[0111] Grid type Total number of cells Number of overlapping cells Percentage of overlapping area Risk point group 200 130 65％ Rock mass deformation boundary 180 130 -

[0112] As shown in Table 5, the grid overlap between the risk point group and the rock mass deformation boundary is counted.

[0113] S502: Call the overlap area ratio to determine whether it exceeds the spatial overlap area ratio critical threshold, if it meets, based on the original boundary curve node coordinates, a cubic spline interpolation function is used to insert new nodes between adjacent nodes, and the first derivative of the interpolated curve is continuous, to generate an expanded boundary curve;

[0114] The spatial overlap area proportion critical threshold is 60% obtained by correlating disaster data statistics and rock force failure probability;

[0115] Determine the overlap area proportion R overlap =0.65 whether it exceeds the critical threshold T=0.6, the comparison operation 0.65>0.6 is true, boundary expansion is performed, and the original vector boundary node coordinate sequence is obtained

[0116] {(x0,y0),(x1,y1),…,(x n ,y n )}={(100,80),(110,85),(120,90),(135,95),(150,85)}, new nodes are inserted between adjacent nodes, for example, two new nodes are inserted between nodes (x1,y1)=(110,85) and (x2,y2)=(120,90), the coordinates of the intermediate points are calculated using a cubic spline interpolation function, the first derivative of the curve is continuous at the nodes after interpolation, and the parameter t is distributed along the curve, and the node spacing is the chord length

[0117] For example, the node spacing The interpolation function S(t)=a(t-t i ) 3 +b(t-t i ) 2 +c(t-t i )+d, where t i is the node parameter, the coefficients a, b, c, and d are determined by the boundary conditions (node coordinates and first derivative continuity), after inserting new nodes, the number of boundary curve nodes increases, for example, the original 5 nodes are expanded to 9 nodes, forming a more detailed boundary description, and an expanded boundary curve is generated.

[0118] Cubic spline interpolation conditions: at node i, 1.S i (t i )=S i+1 (t i )(position continuity) 2.S′ i (t i )=S′ i+1 (t i )(first derivative continuity) 3.S″ i (t i )=S″ i+1 (t i) (second derivative continuous, if using natural spline) select the interpolation between nodes (110, 85) and (120, 90), parameter t in the interval [0, 1], calculate the coordinates at the middle point t = 0.5, assume that through solving the interpolation function to get new nodes (115, 87.5) and (117.5, 88.75).

[0119] Formula: first derivative S'1 at node (110, 85) = (10, 5), first derivative S'2 at node (120, 90) = (15, 5), cubic spline interpolation function coefficients are solved by matrix equation, the coordinates of the interpolation point t = 0.5 are:

[0120] After inserting new nodes, the extended boundary curve is more in line with the actual deformation gradient, providing accurate spatial range for subsequent risk point clustering.

[0121] S503: call the extended boundary curve, extract the coordinate range of the extended area, superimpose the extended area and the risk point group coordinates, perform neighborhood radius coverage detection and minimum point number screening on the risk point group in the extended area based on the density clustering algorithm, merge the point group clusters that meet the conditions, and generate a mine thermal coupling risk map.

[0122] Call the extended boundary curve node coordinates to generate the spatial range of the extended area, for example, the extended boundary coverage coordinate range is x ∈ [95, 155], y ∈ [75, 100], rasterize the extended area, extract all risk point group coordinates within it, for example, risk points (120, 85), (125, 90), (130, 88), (140, 92), apply the density clustering algorithm (DBSCAN), set the neighborhood radius ∈ = 5 meters, the minimum point number N min = 3, traverse each risk point, calculate the number of points in its ∈-neighborhood, for example, the neighborhood of point (125, 90) contains points (120, 85), (130, 88), (125, 90), the point number 3 meets N min , mark it as a core point, merge all density reachable point group clusters, for example, merge {(120, 85), (125, 90), (130, 88)} into a cluster, and generate a mine thermal coupling risk map.

Claims

1. A mine safety production risk prevention and control method based on multi-source data fusion, characterized in that, The method comprises the following steps: S1: acquiring triaxial rock mass displacement rate parameters through a microseismic sensor, and collecting environmental parameters within a corresponding time through a temperature and humidity sensor, constructing a coupling parameter set, calling a Bayesian linear regression model to perform coupling parameter set anomaly detection and extract displacement points respectively, and constructing a rock mass deformation active zone distribution map; S2: determining a displacement area according to the rock mass deformation active zone distribution map, performing difference statistics on the airflow velocity gradient around the displacement area, combining the rock mass displacement rate to extract the rock mass main fracture propagation direction as a boundary condition correction term of a finite volume method model, and constructing a non-steady airflow path curve graph; S3: calculating the curvature change rate of the airflow path curve in the non-steady airflow path curve graph, judging whether the curvature change rates of three consecutive sections exceed 1.5 times of the average value, and positioning a local mutation area of the airflow path according to the judgment result; S4: acquiring the temperature gradient thermal strain rate and axial strain rate of the rock mass material in the local mutation area of the airflow path, generating a corresponding time sequence smooth curve by using a moving average method, identifying a rock mass risk coordinate point group in which the continuous period increments are all positive values in the curve, and constructing a local thermal amplitude trend graph; S5: calling the risk point group of the local thermal amplitude trend graph, calculating the spatial overlap area proportion of the risk point group and a boundary curve of the rock mass deformation active zone distribution map, and if the overlap area proportion is greater than 60%, fitting the boundary curve extension by using a cubic spline interpolation function, and redrawing a new area to generate a mine thermal coupling risk map; The mine thermal coupling risk map comprises a spatial overlap area proportion, a cubic spline interpolation boundary parameter, an extended area coordinate set, and a coupling risk grade partition.

2. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 1, characterized in that: The rock mass deformation active zone distribution map comprises displacement rate anomaly points, environmental parameters, and coupling coefficients, the non-steady airflow path curve graph specifically comprises a velocity gradient difference value, a rock mass main fracture propagation angle, and a corrected boundary condition parameter, the local mutation area of the airflow path specifically refers to a curvature change rate and an airflow velocity anomaly area, and the local thermal amplitude trend graph comprises a temperature strain smooth value, an axial strain increment sequence, and a risk point spatial coordinate.

3. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 1, characterized in that: The specific steps of acquiring triaxial rock mass displacement rate parameters through a microseismic sensor, and collecting environmental parameters within a corresponding time through a temperature and humidity sensor, constructing a coupling parameter set, calling a Bayesian linear regression model to perform coupling parameter set anomaly detection and extract displacement points, and constructing a rock mass deformation active zone distribution map are as follows: S101: acquiring triaxial displacement rate data of a microseismic sensor and airflow velocity, temperature gradient, and humidity saturation parameters collected by a temperature and humidity sensor, and processing different dimension parameters by using a minimum-maximum normalization method, performing spatial coordinate matching on the four types of parameters based on a unified time stamp, filling in missing values by using a linear interpolation method, and generating a coupling parameter set; S102: acquiring displacement rate, airflow velocity, temperature gradient, and humidity saturation parameters of the coupling parameter set, constructing a linear regression equation through a Bayesian linear regression model, calculating residual absolute values and a 95% confidence interval, screening measurement point coordinates whose residuals exceed the upper limit of the confidence interval, and generating an anomaly detection result set; The Bayesian linear regression model sets the prior distribution of regression coefficients and noise, combines observation data to deduce the posterior distribution, generates the predicted value and confidence interval of the displacement rate, and is used for residual analysis and abnormal point screening; S103: calling the abnormal detection result set, extracting the coordinates of the measurement points with a displacement rate of greater than or equal to 0.5 mm / h, using a spatial neighborhood clustering algorithm to perform density clustering on the measurement points, merging adjacent points with a spacing less than a neighborhood radius, determining the neighborhood radius by using a K nearest neighbor method according to the average distance from a measurement point to a Kth nearest neighbor point, and generating a rock mass deformation active zone distribution map.

4. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 3, characterized in that: According to the rock mass deformation active zone distribution map, the displacement region is determined, the airflow velocity gradient around the displacement region is differentially statistically analyzed, the rock mass main fracture propagation direction is extracted as a boundary condition correction term of a finite volume method model in combination with the rock mass displacement rate, and the specific steps of constructing a non-steady airflow path curve diagram are as follows: S201: calling the rock mass deformation active zone distribution map, extracting the coordinates of the measurement points around the displacement region, calculating the airflow velocity gradient based on the Euclidean distance between the measurement points and the airflow velocity difference between adjacent measurement points, and generating a gradient statistical result; S202: calling the gradient statistical result, comparing the airflow velocity gradient of each measurement point with a critical threshold value of airflow velocity gradient perturbation error item by item, screening the measurement point coordinates with an airflow velocity gradient exceeding the critical threshold value of airflow velocity gradient perturbation error, and generating an over-limit region coordinate set; The critical threshold value of the airflow velocity gradient perturbation error is determined by parameter sensitivity analysis in the finite volume method, when the gradient exceeds 0.3, and the perturbation error of the airflow to the rock mass displacement direction exceeds 5%, the target is taken as a trigger point for boundary condition correction; The airflow velocity gradient difference statistics are determined by calculating the Euclidean distance mean value of the first-order derivative of the airflow velocity vector in the neighborhood; S203: calling the over-limit region coordinate set, calculating the feature vector of the rock mass displacement rate by using a geometric direction statistical method, extracting the feature vector direction corresponding to the maximum eigenvalue, taking the rock mass main fracture propagation direction as the normal correction amount of the boundary condition of the finite volume method, and generating a non-steady airflow path curve diagram; The correction term of the finite volume method is the projection component of the original normal velocity component along the rock mass main fracture propagation direction.

5. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 4, characterized in that: The specific steps of calculating the curvature change rate of the airflow path curve in the non-steady airflow path curve diagram and judging whether the curvature change rates of three consecutive segments exceed 1.5 times the mean value according to the judgment result to locate the local mutation region of the airflow path are as follows: S301: obtaining a discrete coordinate point set of the non-steady airflow path curve diagram, parameterizing the curve by using a cubic spline interpolation method, the condition of the cubic spline interpolation being that the first-order derivative and the second-order derivative of the curve at each interpolation node are continuous, calculating the curvature value of the non-steady airflow path curve diagram point by point by using a parameterized cubic spline equation, generating a curvature value sequence, and performing a difference operation on the curvatures of adjacent sampling points to generate a curvature change rate sequence; S302: Call the curvature change rate sequence, calculate the mean of all data points in the curvature change rate sequence, compare the values of the three adjacent data points in the sequence with the mean, and if the absolute values of the three data points are all more than 1.5 times the mean, mark the interval and generate a super-mean marked sequence; The curvature change rate is calculated by using a first-order difference method, that is, the curvature values of adjacent sampling points are differentiated to approximate the curvature change rate by the difference between the curvature of the current sampling point and the curvature of the previous sampling point; S303: Based on the super-mean marked sequence, extract the starting and ending coordinate point indexes of the marked interval, merge the intervals with adjacent intervals less than five sampling points, determine the boundary range of the curvature mutation region in the airflow path, and locate the local mutation region of the airflow path.

6. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 5, characterized in that: Parametric planar curves At parameter The processing is performed using the formula: ； Computing curvature ; in, Parametric curves The parameters, Is the curve in the parameter The x-coordinate function at that location Is the curve in the parameter The y-coordinate function at the location, parameters and They are curves and Coordinates with respect to parameters The first derivative represents the tangent vector component of the curve at that point, and the parameter... and It is the second derivative, related to the acceleration or bending pattern of the curve, in the molecule. It is a two-dimensional form of the cross product of the tangent vector and the acceleration vector, with the denominator being... It is the speed raised to the power of 3 / 2.

7. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 6, characterized in that: The specific steps of obtaining the temperature gradient thermal strain rate and axial strain rate of the rock mass material in the local mutation region of the airflow path, generating the corresponding time sequence smooth curve by using the moving average method, identifying the rock mass risk coordinate point group in which the continuous period increments are all positive values, and constructing the local thermal force amplitude trend chart are as follows: S401: Call the local mutation region coordinates of the airflow path, extract the temperature gradient thermal strain rate and axial strain rate original sampling point data of the rock mass material based on the coordinate range, align the time sequence data of the two kinds of strain rates by interpolation, map the strain rate values to the 0-1 interval by using linear normalization method, and generate time-aligned strain rate data set; S402: According to the time-aligned strain rate data set, a sliding window is established for the temperature gradient thermal strain rate and the axial strain rate respectively, the window covers the front two and the back two data points, the arithmetic mean of the five data points in the window is calculated to replace the original center point value, and the smoothed strain rate curve set is generated; The window size definition of the moving average method is based on the signal main period setting reference, combined with the noise standard deviation evaluation to set the lower limit, and verified by autocorrelation to ensure trend continuity, and dynamically adjusted to balance noise suppression and trend capture; S403: Based on each time node of the smoothed strain rate curve set, the difference value of the temperature gradient thermal strain rate and the axial strain rate is calculated, the coordinate points with a difference value of more than zero for three consecutive periods are extracted, and the coordinate points that meet the condition are projected to a two-dimensional grid with time axis horizontally and strain rate amplitude vertically according to the spatial coordinates, and a local thermal force amplitude trend chart is constructed.

8. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 7, characterized in that: Call the risk point group of the local thermal force amplitude trend chart, calculate the spatial overlap area ratio of the boundary curve of the rock mass deformation active zone distribution map, if the overlap area ratio is greater than 60%, use the cubic spline interpolation function to fit the boundary curve expansion, and the specific steps of redrawing the mine thermal coupling risk map for the new area are as follows: S501: Obtain the risk point group of the local thermal force amplitude trend chart and the boundary curve of the rock mass deformation active zone distribution map, convert the vector boundary to raster data, count the number of pixels of the intersection of the raster where the risk point group is located and the boundary raster, calculate the proportion of the intersection pixel number to the total risk point group pixel number, and generate the overlap area ratio; S502: Call the overlap area ratio, judge whether it exceeds the space overlap area ratio critical threshold, if it meets, based on the original boundary curve node coordinates, use cubic spline interpolation function to insert new nodes between adjacent nodes, constrain the first derivative of the interpolated curve to be continuous, and generate an expanded boundary curve; The space overlap area ratio critical threshold is 60% of the space overlap area critical threshold obtained by associating the disaster data statistics with the rock force failure probability; S503: Call the expanded boundary curve, extract the extended area coordinate range, superimpose the extended area and the risk point group coordinates, perform neighborhood radius coverage detection and minimum point number screening on the risk point group in the extended area based on the density clustering algorithm, merge the point group clusters that meet the conditions, and generate a mine thermal coupling risk map.

9. The mine safety production risk prevention and control method based on multi-source data fusion according to claim 8, characterized in that: The proportion of the overlapping area of the risk point group and the deformation boundary is calculated by a logistic regression model, and the formula is: ; and the rock force failure probability is calculated ; wherein, is the ratio of the overlapping area of the risk point and the boundary, representing the statistical result of the disaster data, is the logarithmic odds of the baseline failure probability, represents the influence coefficient of the overlapping area on the failure probability, reflecting the increase amplitude of the logarithmic odds of the failure probability when the overlapping area increases by 1%.

Images