Intelligent prediction method for space-time evolution of coal mining subsidence

Abstract

The invention relates to the technical field of coal mining safety monitoring, and discloses a coal mining subsidence space-time evolution intelligent prediction method, which comprises the following steps: acquiring multi-dimensional geological dynamic data such as surface displacement, underground water pressure and rock stratum stress of a coal mining area, extracting a space-time characteristic sequence by using a GCLSTM algorithm, and calculating a space-time fluctuation coefficient; extracting a principal component vector through spectral decomposition, mapping to obtain a space-time correlation weight, dynamically correcting an original prediction value to generate a fusion prediction value, and judging a subsidence space-time evolution trend; constructing a three-dimensional subsidence curved surface model to identify a high-risk area and generate graded early warning; and establishing a joint optimization model to solve an optimal mining path, and generating a subsidence risk minimization scheme. According to the method, intelligent analysis of multi-dimensional data fusion and spatial-temporal characteristics is realized, the prediction accuracy is improved, an effective means is provided for coal mining subsidence prevention and control and safe mining, and the method has a remarkable engineering application value.

Description

Technical Field

[0001] The present invention relates to the technical field of coal mining safety monitoring, and in particular to an intelligent prediction method for the spatiotemporal evolution of coal mining subsidence. Background Art

[0002] Mining subsidence, a common geological hazard in coal mining, poses a serious threat to surface structures, the ecological environment, and groundwater resources. With the increasing depth and intensity of coal mining, the complexity and uncertainty of subsidence hazards have increased significantly. Traditional subsidence prediction methods, unable to effectively capture the spatiotemporal coupling characteristics of multi-dimensional geological data, lack accuracy and are unable to meet the requirements of safe and efficient modern coal mining.

[0003] From a technological perspective, early subsidence predictions primarily relied on empirical formulas and single-sensor monitoring data, such as surface movement prediction models based on the probability integral method. These methods can only describe subsidence patterns under simple geological conditions and are unable to adapt to the demands of complex geological structures and dynamic mining processes. With the advancement of sensor technology, the fusion of multi-source monitoring data has become a trend. However, traditional data fusion methods (such as weighted averaging and Kalman filtering) are mostly based on linear assumptions and have difficulty capturing the nonlinear spatiotemporal correlations between geological data. For example, multidimensional data such as surface displacement, groundwater pressure, and rock stress interact in complex temporal and spatial ways, and linear models cannot accurately capture their dynamic coupling relationships, resulting in significant deviations between predicted results and the actual subsidence evolution process.

[0004] After deep learning technology was introduced into the field of geological disaster prediction, recurrent neural networks (RNN) and its variants (such as LSTM and GRU) have been widely used due to their ability to model time series features. However, traditional LSTM models have limitations in processing spatial features and cannot effectively capture the local correlation of geological data on spatial grids. Although gated convolutional networks (Gated Convolutional Networks) have shown advantages in extracting spatiotemporal features from images and videos, when applied alone to geological data prediction, it is difficult to fully exploit the long-term dependency characteristics of time series. In addition, when dealing with the weight allocation of multi-dimensional data, existing methods usually use fixed weights or simple empirical assignments, ignoring the dynamic importance differences of different geological factors in different spatiotemporal scenarios, resulting in insufficient generalization ability of the prediction model.

[0005] When it comes to subsidence risk assessment and mining plan optimization, traditional methods often analyze prediction results and mining schedules independently, lacking a dynamic coupling mechanism. For example, mining path planning is typically based on static geological models, failing to consider the feedback adjustment of mining sequences based on real-time subsidence prediction results, making it difficult to achieve active control of subsidence risk. Furthermore, existing three-dimensional subsidence model construction methods (such as polynomial interpolation and triangulated network modeling) are prone to surface distortion when processing high-noise discrete data, affecting the accuracy of identifying subsidence risk areas. Summary of the Invention

[0006] The purpose of the present invention is to provide an intelligent prediction method for the spatiotemporal evolution of coal mining subsidence to solve the problems raised in the above background technology.

[0007] To achieve the above objectives, the present invention provides the following technical solution: an intelligent prediction method for the spatiotemporal evolution of coal mining subsidence, the method comprising:

[0008] Collect multi-dimensional geological dynamic data of the coal mining area, wherein the geological dynamic data includes surface displacement, groundwater pressure and rock formation stress;

[0009] The GCLSTM algorithm is used to extract spatiotemporal features from geological dynamic data to obtain spatiotemporal feature sequences, and the spatiotemporal fluctuation coefficient is calculated based on the gradient changes of the spatiotemporal feature sequences.

[0010] The covariance matrix of the spatiotemporal fluctuation coefficients of different dimensions is spectrally decomposed to extract the principal component vectors, which are then nonlinearly mapped using the hyperbolic tangent function to obtain the spatiotemporal correlation weights of each dimension.

[0011] The original prediction value of the spatiotemporal feature sequence is dynamically modified according to the spatiotemporal correlation weight to generate a fused prediction value, and the spatiotemporal evolution trend of coal mining subsidence is determined based on the fused prediction value;

[0012] The fusion prediction value F′ satisfies the relationship:

[0013] F′=F×tanh(SC×W)

[0014] Where F represents the original prediction value, SC represents the spatiotemporal fluctuation coefficient, W represents the spatiotemporal correlation weight, and tanh() represents the hyperbolic tangent function.

[0015] Preferably, the spatiotemporal feature extraction includes:

[0016] Intercept geological dynamic data in a preset time window, construct a space-time tensor and perform normalization processing;

[0017] A gated convolutional layer is used to capture local spatiotemporal features of the spatiotemporal tensor and output the Hadamard product of the gated weights and the convolutional features;

[0018] The gated convolution result is input into the bidirectional LSTM unit to extract the forward and backward long-term dependency features respectively, and the feature sequence is weightedly fused through the attention mechanism to generate a spatiotemporal feature sequence.

[0019] Preferably, the calculation of the spatiotemporal fluctuation coefficient includes:

[0020] Select the spatiotemporal feature sequence at any moment as the benchmark feature, and calculate the gradient modulus of the feature vectors at the next n consecutive moments;

[0021] The ratio of the mean value of the gradient modulus to the Euclidean distance of the benchmark feature is used as the time dimension fluctuation;

[0022] Calculate the KL divergence between the benchmark feature and the adjacent spatial grid feature, and take its reciprocal as the spatial dimension fluctuation;

[0023] The geometric mean of the time dimension fluctuation and the space dimension fluctuation is taken as the space-time fluctuation coefficient.

[0024] Preferably, the principal component vector extraction includes:

[0025] Construct the covariance matrix of spatiotemporal fluctuation coefficients of different dimensions, perform singular value decomposition on it and obtain the eigenvalue diagonal matrix;

[0026] Select eigenvectors whose eigenvalue cumulative contribution rate exceeds a preset threshold to form the principal component space;

[0027] Project the covariance matrix into the principal component space to obtain the principal component vector after dimensionality reduction.

[0028] Preferably, the spatiotemporal association weight calculation includes:

[0029] Normalize the principal component vector and calculate its cosine similarity with the preset reference vector;

[0030] The cosine similarity is input into the multi-layer perceptron, and the initial weight is generated after mapping through the fully connected layer and activation function;

[0031] The initial weights are temporally smoothed through the gated recurrent unit to output the spatiotemporal correlation weights of each dimension.

[0032] Preferably, the generation of the fusion prediction value includes:

[0033] Perform tensor multiplication of the spatiotemporal correlation weight and the original prediction value to obtain the weight-corrected prediction value;

[0034] Calculate the residual between the weighted corrected prediction value and the original prediction value, and dynamically compensate the residual through the Kalman filter;

[0035] The compensated residual is added to the original prediction value to generate a fused prediction value.

[0036] Preferably, the method further comprises:

[0037] Based on the fusion prediction values, a three-dimensional subsidence surface model is constructed to extract the curvature extreme points and gradient mutation areas of the surface;

[0038] When the number of extreme curvature points exceeds the preset threshold or the area proportion of the gradient mutation region is greater than the critical value, it is judged as a high-risk area for subsidence and a graded warning signal is generated.

[0039] Preferably, the three-dimensional subsidence surface model construction includes:

[0040] The fused prediction value is mapped to a three-dimensional coordinate system according to the spatial grid coordinates to generate discrete point cloud data;

[0041] The radial basis function interpolation method is used to fit the surface of discrete point clouds, and the Gaussian curvature and mean curvature of the fitting surface are calculated;

[0042] The curvature distribution is smoothed by moving least squares method to eliminate noise interference.

[0043] Preferably, the method further comprises:

[0044] A joint optimization model of geological dynamic data and mining progress parameters was established, and the alternating direction multiplier method was used to solve the optimal mining path;

[0045] The optimal mining path is dynamically coupled with the fused prediction value to generate a mining plan that minimizes the subsidence risk.

[0046] Preferably, the joint optimization model solution includes:

[0047] The objective function is defined as the weighted square error sum of the subsidence prediction value and the mining progress, and the constraint condition is the rock formation stress safety threshold;

[0048] Perform Lagrangian dual transformation on the objective function and decompose it into mining progress subproblem and subsidence prediction subproblem;

[0049] The two sub-problems are solved alternately and iteratively until convergence, and the optimal mining path that meets the constraints is output.

[0050] Compared with the prior art, the present invention has the following beneficial effects:

[0051] In terms of data processing, the collection of multi-dimensional geological dynamic data, including surface displacement, groundwater pressure, and rock layer stress, overcomes the limitations of traditional single-dimensional data. Multi-dimensional data can comprehensively reflect the geological state of the mining area, providing a rich information foundation for accurate prediction. The GCLSTM algorithm extracts spatiotemporal features from the geological dynamic data. Combining gated convolutional layers and bidirectional LSTM units, it effectively captures local spatiotemporal features and long-term dependencies. At the same time, the attention mechanism is used to weightedly fuse feature sequences, making the extracted spatiotemporal features more accurate and comprehensive, better reflecting the spatiotemporal evolution of geological data.

[0052] In terms of feature analysis, by calculating the spatiotemporal fluctuation coefficient and extracting the principal component vectors, we deeply explore the spatiotemporal correlation characteristics between data of different dimensions. The spatiotemporal fluctuation coefficient comprehensively considers the fluctuations in both the temporal and spatial dimensions, accurately reflecting the degree of dynamic change in the spatiotemporal feature sequence. Spectral decomposition of the covariance matrix and extraction of the principal component vectors achieve dimensionality reduction for high-dimensional data, preserving key information while reducing computational complexity. Spatiotemporal correlation weights are obtained by nonlinearly mapping the principal component vectors using the hyperbolic tangent function. This accurately describes the nonlinear correlation between data of different dimensions, providing a scientific basis for subsequent corrections to the predicted values.

[0053] In predictive model construction, the original prediction values ​​of the spatiotemporal feature sequence are dynamically modified based on spatiotemporal correlation weights to generate fused prediction values. Dynamic compensation for residual errors is achieved through tensor multiplication and a Kalman filter, effectively improving the accuracy and reliability of the prediction results. The fused prediction values ​​more realistically reflect the actual evolution of coal mining subsidence, providing a reliable basis for determining subsidence risk.

[0054] In terms of risk warning and mining optimization, a three-dimensional subsidence surface model is constructed based on the fused prediction values. By extracting curvature extreme points and gradient mutation areas, it can accurately identify high-risk subsidence areas and generate graded warning signals, achieving real-time warning and risk control for coal mining subsidence. At the same time, a joint optimization model of geological dynamic data and mining progress parameters is established. The alternating direction multiplier method is used to solve the optimal mining path, which is dynamically coupled with the fused prediction values ​​to generate a mining plan that minimizes subsidence risk. This organically combines prediction results with mining path planning, effectively reducing coal mining subsidence risks and improving the safety and economic efficiency of coal resource mining. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 This is a working principle diagram of the intelligent prediction method for the spatiotemporal evolution of coal mining subsidence according to the present invention;

[0056] Figure 2 Design diagram for spatiotemporal feature extraction;

[0057] Figure 3 This is the design drawing of the subsidence early warning system. DETAILED DESCRIPTION

[0058] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0059] See also Figure 1-Figure 3 The present invention relates to an intelligent prediction method for the spatiotemporal evolution of coal mining subsidence, and the specific implementation steps are as follows:

[0060] Multi-dimensional geological dynamic data collection is conducted in the coal mining area, covering three core data types: surface displacement, groundwater pressure, and rock formation stress. Surface displacement is monitored in real time using the Beidou satellite differential positioning system, with a sampling interval of 10 minutes and centimeter-level positioning accuracy. Groundwater pressure is monitored using a distributed fiber-optic sensing network, with sensor nodes deployed in the aquifers, acquiring water pressure data every hour. Rock formation stress is collected using strain gauges embedded in the coal seam roof, with a sampling frequency of once per minute. All data is aggregated in real time to a data center via a wireless transmission module, forming a multi-dimensional dynamic dataset.

[0061] The GCLSTM algorithm is used to extract spatiotemporal features from the collected geological dynamic data. Specifically, data is captured within a preset 24-hour time window. The multidimensional data within each time window is constructed into a three-dimensional spatiotemporal tensor (time dimension × spatial grid dimension × feature dimension). The data is then normalized to the range [-1, 1]. Subsequently, a gated convolutional layer is used to capture local spatiotemporal features from the spatiotemporal tensor. This layer contains 32 3×3 convolution kernels. A gating mechanism generates a gating weight vector, which is then Hadamard-producted with the convolutional features to output a local feature response. The gated convolution result is input into a bidirectional LSTM unit, which contains 128 hidden neurons. This unit extracts long-term dependency features in the forward and backward directions, respectively. The forward and backward feature sequences are weighted and fused using an attention mechanism to generate a feature sequence containing spatiotemporal coupling information.

[0062] The spatiotemporal feature sequence at the initial moment is selected as the benchmark feature, and the gradient modulus of the feature vectors at the subsequent n consecutive moments (n=10) is calculated. The gradient calculation adopts the central difference method. The ratio of the mean of the gradient modulus to the Euclidean distance of the benchmark feature is defined as the time dimension fluctuation, which reflects the rate of change of the time series. At the same time, the KL divergence between the benchmark feature and the adjacent spatial grid feature is calculated, and its reciprocal is taken as the spatial dimension fluctuation to measure the feature difference between adjacent spatial regions. Finally, the spatiotemporal fluctuation coefficient is determined by the geometric mean of the time dimension fluctuation and the spatial dimension fluctuation. The calculation formula is:

[0063]

[0064] The covariance matrix of the spatiotemporal fluctuation coefficients of different dimensions is constructed, and the singular value decomposition is performed on it to obtain the eigenvalue diagonal matrix and the eigenvector matrix. The eigenvectors with a cumulative contribution rate of more than 85% of the eigenvalues ​​are selected to form the principal component space, and the covariance matrix is ​​projected onto the principal component space to obtain the principal component vector after dimensionality reduction. After the principal component vector is standardized, its cosine similarity with the preset reference vector (a typical eigenvector obtained from historical sinking data statistics) is calculated, and the similarity value is input into a multilayer perceptron with two hidden layers (64 neurons in each layer, with ReLU activation function). After mapping through the fully connected layer and the hyperbolic tangent activation function, the initial weights are generated. Finally, the initial weights are time-series smoothed by the gated recurrent unit (GRU). The GRU contains 32 hidden units and outputs the spatiotemporal correlation weights of each dimension.

[0065] The original prediction value of the spatiotemporal feature sequence is dynamically corrected based on the spatiotemporal correlation weight. First, the spatiotemporal correlation weight is tensor-multiplied with the original prediction value to obtain the weight-corrected prediction value; then, the residual between the weight-corrected prediction value and the original prediction value is calculated, and the residual is dynamically compensated through the Kalman filter. The state transfer matrix and observation matrix of the Kalman filter are determined based on historical data training; finally, the compensated residual is superimposed on the original prediction value to generate the fused prediction value F', which is calculated as follows:

[0066] F′=F×tanh(SC×W)

[0067] Based on the fused prediction values, the temporal evolution trend (such as the change in subsidence rate) and spatial distribution trend (such as the migration of subsidence centers) of coal mining subsidence are determined through time series trend analysis algorithms (such as the Holt-Winters model) and spatial interpolation algorithms (such as Kriging interpolation).

[0068] The present invention will be further described below in conjunction with Examples 1 to 5:

[0069] Example 1:

[0070] Based on the above overall scheme, the spatiotemporal feature extraction process further includes the following specific implementation methods:

[0071] Geological dynamic data is captured within a preset time window. The duration of the preset time window can be determined based on the geological conditions and mining progress of the coal mining area. For example, it can be set to 24 hours, with the geological dynamic data within 24 consecutive hours being used as an analysis unit. Within this time window, three types of data are synchronously collected: surface displacement, groundwater pressure, and rock formation stress. Surface displacement data is acquired via the Beidou satellite differential positioning system at a sampling interval of 10 minutes, with each time point corresponding to the displacement coordinates of each spatial grid within the coal mining area. Groundwater pressure data is collected via distributed fiber optic sensing nodes deployed in the aquifer at a sampling interval of 1 hour, with each node recording the water pressure value at its location. Rock formation stress data is acquired via strain gauges embedded in the coal seam roof at a sampling interval of 1 minute, with each strain gauge monitoring stress changes in the rock formation in real time.

[0072] The collected data must be constructed as a spatiotemporal tensor. The dimensions of a spatiotemporal tensor are defined as (T × N × D), where T is the time step, corresponding to the number of sampling points within a preset time window. For example, a 24-hour period contains 144 10-minute time points (surface displacement data), 24 1-hour time points (groundwater pressure data), or 1440 1-minute time points (rock stress data). In practical applications, the time steps of data of different dimensions must be unified through interpolation or downsampling. For example, if all are converted to a 10-minute step, then T = 144. N is the number of spatial grids. The coal mining area can be divided into a regular square grid based on geographic coordinates, such as a 100 × 100 grid, in which case N = 10,000. D is the feature dimension, corresponding to surface displacement (displacement components in the x, y, and z directions), groundwater pressure (water pressure value), and rock stress (vertical and horizontal stress components), totaling six feature dimensions (D = 6). Therefore, the data in each time window constitute a three-dimensional space-time tensor of 144×10000×6, representing the multi-dimensional geological characteristics at different time points and different spatial grids.

[0073] After constructing the spatiotemporal tensor, it needs to be normalized to eliminate the dimensional differences of different feature dimensions. Normalization uses the Z-score normalization method, and the specific formula is:

[0074]

[0075] Where X is the original eigenvalue, μ is the mean of the feature dimension within the time window, and σ is the standard deviation. This process normalizes the data of each dimension to a standard normal distribution with a mean of 0 and a standard deviation of 1, facilitating feature learning in subsequent neural networks.

[0076] A gated convolutional layer is used to capture local spatiotemporal features of the normalized spatiotemporal tensor. The gated convolutional layer contains 32 3×3 convolution kernels. The convolution operation is performed on the spatial dimension (N dimensions) with a stride of 1. Same padding is used to ensure that the size of the convolved feature map matches the input. Each convolution kernel generates a two-dimensional feature map (T×N) representing the characteristic response of the corresponding local spatial region. The gating mechanism is implemented using an independent gating weight vector, generated using a sigmoid activation function and ranging from [0, 1]. This vector is then multiplied (element-wise) with the convolution feature map to suppress invalid or noisy information in the convolution features while retaining key local spatiotemporal features. For example, for abnormal fluctuations in surface displacement characteristics caused by sensor errors, the gating weights can automatically learn to lower values, reducing the impact of these features.

[0077] The output of the gated convolutional layer is input to a bidirectional LSTM unit. The bidirectional LSTM consists of a forward LSTM and a backward LSTM, each containing 128 hidden neurons. The forward LSTM processes the feature sequence in chronological order (t = 1 to t = T), capturing long-term dependencies in the past; the backward LSTM processes the feature sequence in reverse chronological order (t = T to t = 1), capturing dependencies in the future. The hidden states of the two are fused through a concatenation operation to form a feature vector sequence containing bidirectional temporal information. For example, for a rock formation stress sequence in a spatial grid, the forward LSTM can learn the cumulative impact of mining activities on rock formation stress, while the backward LSTM can learn the potential impact of subsequent mining plans on the current stress state. This bidirectional fusion more comprehensively reflects the temporal dependence of stress changes.

[0078] After the bidirectional LSTM outputs the feature sequence, the feature sequence is weighted and fused through the attention mechanism. The attention mechanism adopts the additive attention model. The specific calculation process is as follows: First, the hidden state sequence h of the bidirectional LSTM is t (t=1,2,...,T) is linearly transformed with the trainable query vector q, and the energy value e is generated through the tanh activation function t :

[0079]

[0080] Among them, W h and W q is the weight matrix, w h is the output weight vector, which is the parameter learned during the model training process. Then, the energy value is normalized by the softmax function to obtain the attention weight α of each time step t :

[0081]

[0082] The attention weight reflects the importance of each time-step feature to the final spatiotemporal feature sequence. For example, during periods of intense mining activity, the corresponding attention weight increases significantly, indicating that the geological characteristics of that period have a higher contribution to subsidence prediction. Finally, the attention weight is weighted and summed with the hidden state sequence to generate the final spatiotemporal feature sequence s:

[0083]

[0084] This spatiotemporal feature sequence combines the long-term dependence characteristics of the time dimension and the local correlation characteristics of the space dimension, and can effectively characterize the spatiotemporal coupling characteristics of coal mining subsidence.

[0085] In the implementation, the parameters of the gated convolutional layers and bidirectional LSTM units were trained using a backpropagation algorithm. The mean squared error (MSE) loss function was used, and the Adam optimizer was selected as the optimizer. The learning rate was initialized to 0.001 and adjusted using an exponential decay strategy. During model training, historical coal mining subsidence data was used as the training set, and cross-validation was used to divide the training and validation sets to ensure the model's generalization ability.

[0086] Example 2:

[0087] The specific implementation of the calculation of the spatiotemporal fluctuation coefficient is as follows:

[0088] In the calculation of the time dimension fluctuation, the spatiotemporal feature sequence at any moment needs to be selected as the benchmark feature. Usually, the feature sequence f0 at the initial moment of the time window (such as t = 0) is selected as the benchmark. This feature sequence is the vector corresponding to the initial moment in the spatiotemporal feature sequence output by the spatiotemporal feature extraction step, and the dimension is consistent with the feature dimension (such as a fusion vector containing multi-dimensional features such as surface displacement, groundwater pressure, and rock formation stress). Subsequently, the gradient modulus of the feature vectors at n consecutive moments after the benchmark feature is calculated, where n is a preset time step parameter, which can be determined according to the time evolution scale of coal mining subsidence. For example, n = 10 is set (i.e., analyzing the feature changes in the subsequent 10 time steps).

[0089] The gradient modulus is calculated using the L2 norm (Euclidean distance), and the specific formula is:

[0090]

[0091] Among them, f t represents the eigenvector at time t, f t，iis the value of the i-th feature dimension of the vector (i = 1, 2, ..., D, where D is the total number of feature dimensions). This calculation reflects the absolute change in the feature vectors at adjacent time steps. The gradient modulus at n consecutive moments can be used to characterize the fluctuation amplitude of the feature sequence in the time dimension.

[0092] After calculating the gradient modulus of n consecutive moments, take its mean as the average rate of change in the time dimension. The formula is:

[0093]

[0094] In order to eliminate the influence of the amplitude of the benchmark feature itself, the average gradient modulus and the Euclidean distance of the benchmark feature are normalized to obtain the time dimension fluctuation amount, which is:

[0095]

[0096] This fluctuation quantity characterizes the relative change intensity of the characteristic sequence in the time dimension. The value range is non-negative. The larger the value, the more intense the fluctuation of the time series.

[0097] In the calculation of spatial fluctuations, the KL divergence (Kullback-Leibler Divergence) is used to measure the distribution difference between the benchmark feature and the adjacent spatial grid features. The definition of adjacent spatial grids can be determined based on the actual spatial resolution. For example, for a two-dimensional spatial grid partition (row i, column j), the benchmark feature is located at grid (i, j). Its adjacent grids can include four direct neighboring grids (i±1, j) and (i, j±1) above, below, and left, or right, or can be expanded to eight neighboring grids.

[0098] For each adjacent grid (k, l), calculate its eigenvector f k，l The KL divergence with the benchmark feature f0 is:

[0099]

[0100] Among them, f k，l，i is the value of the i-th dimension of the feature vector of the adjacent grid (k, l). KL divergence is an asymmetric metric used to measure the difference between two probability distributions. Here, the values ​​of each dimension of the feature vector are considered to be normalized probability distributions (before actual calculation, the feature vector can be normalized using the softmax function to ensure that the values ​​of each dimension are non-negative and sum to 1). Therefore, a larger KL divergence value indicates a more significant spatial difference between the baseline feature and the adjacent grid feature.

[0101] In order to make the spatial dimension fluctuation amount and the time dimension fluctuation amount have consistent monotonicity (the larger the value, the stronger the fluctuation), the inverse of the KL divergence is taken as the spatial dimension fluctuation amount, and the formula is:

[0102]

[0103] Among them, ∈ is the minimum value (such as 10 -8 ), which is used to avoid the situation where the denominator is zero. For multiple adjacent grids (such as four or eight neighborhoods), the spatial dimension fluctuation can be the average value of the fluctuation of each neighborhood, and the formula is:

[0104]

[0105] Where M is the number of adjacent grids. This fluctuation range is positive. Larger values ​​indicate smaller differences between the baseline characteristics and those of adjacent spatial grids, indicating a more uniform spatial distribution. Conversely, smaller values ​​indicate greater spatial differences, suggesting the presence of areas with potential spatial mutations at risk of subsidence.

[0106] The space-time fluctuation coefficient is determined by the geometric mean of the time dimension fluctuation and the space dimension fluctuation, and the formula is:

[0107]

[0108] The geometric mean calculation method balances the influence of fluctuations in both time and space, preventing extreme values ​​in a single dimension from overly dominating the overall fluctuation coefficient. For example, when the fluctuation in the time dimension is large but the fluctuation in the space dimension is small, the spatiotemporal fluctuation coefficient will not be significantly high, and vice versa. This method more comprehensively reflects the spatiotemporal coupled fluctuation characteristics of geological characteristics.

[0109] In a specific implementation, the value of n in the time dimension and the range of the spatial neighborhood can be used as model hyperparameters, with optimal values ​​determined through sensitivity analysis of historical data. For example, in areas with rapid mining rates, the value of n can be reduced to capture short-term, drastic changes; in areas with complex geological conditions, the spatial neighborhood can be expanded to eight neighborhoods to more comprehensively assess spatial differences. Furthermore, the normalization of feature vectors (such as the softmax used for probability conversion before KL divergence calculation) and the setting of the minimum value ∈ must strictly follow numerical computation specifications in the code implementation to avoid computational overflow or loss of precision.

[0110] Example 3:

[0111] The principal component vector extraction process includes the following specific implementation methods:

[0112] Construct the covariance matrix of spatiotemporal fluctuation coefficients of different dimensions. The spatiotemporal fluctuation coefficient involves multiple dimensions such as surface displacement, groundwater pressure, rock stress, etc. (assuming D dimensions, D ≥ 2). The spatiotemporal fluctuation coefficient of each dimension constitutes a time series. For example, for the surface displacement dimension, each time step corresponds to a fluctuation coefficient value, and the same is true for other dimensions. Assume that there are T time steps in total, and the fluctuation coefficient sequence of each dimension is {sc1(t), sc2(t), ..., sc D (t)}(t=1,2,…,T), then the covariance matrix C is a D×D matrix, where the element C i，j It represents the covariance of the spatiotemporal fluctuation coefficients of the i-th dimension and the j-th dimension, and the calculation formula is:

[0113]

[0114] in, and are the sample means of the volatility coefficients of the i-th and j-th dimensions, respectively. The covariance matrix is ​​used to characterize the linear correlation between the volatility coefficients of different dimensions. The diagonal elements are the variances of each dimension, and the off-diagonal elements are the covariances between any two dimensions.

[0115] Next, perform singular value decomposition (SVD) on the covariance matrix C to obtain three matrices: the left singular matrix U, the singular value diagonal matrix Σ, and the right singular matrix V, namely:

[0116] C=U∑V T

[0117] Among them, U and V are both orthogonal matrices (U T U=1,V T V = I, I is the unit matrix), Σ is a diagonal matrix, the diagonal elements are the singular values ​​of the covariance matrix (i.e., the square roots of the eigenvalues), and are arranged in descending order as σ1≥σ2≥…≥σ D ≥0. Each singular value corresponds to an eigenvalue Characterizes the variance contribution of the corresponding principal component.

[0118] The principal component space is constructed by selecting eigenvectors whose eigenvalue cumulative contribution rate exceeds a preset threshold. The preset threshold is usually set to 80% to 95%. For example, 85% means retaining the principal components that can explain 85% of the variance of the original data. The formula for calculating the cumulative contribution rate is:

[0119]

[0120] Where k is the number of principal components selected (k≤D), and the minimum k value that satisfies the cumulative contribution rate ≥ the preset threshold. The corresponding first k eigenvectors (i.e., the first k column vectors u1, u2, ..., u in the left singular matrix U) k) constitute the basis vectors of the principal component space, each basis vector represents a principal component and reflects the direction of the linear combination of the original dimensions.

[0121] After that, the covariance matrix is ​​projected into the principal component space to obtain the principal component vector after dimensionality reduction. For the original D-dimensional spatiotemporal fluctuation coefficient vector sc=[sc1,sc2,...,sc D ] T , its projection in the principal component space is:

[0122]

[0123] Among them, U k is a D×k matrix composed of the first k columns of the left singular matrix U, pc is the k-dimensional principal component vector, and each element pc i (i=1,2,…,k) represents the score of the i-th principal component, reflecting the projection strength of the original fluctuation coefficient in the direction of the principal component.

[0124] In practice, the construction of the covariance matrix requires careful data centering. Specifically, the mean of the coefficient of fluctuation for each dimension must be subtracted to ensure the accuracy of the covariance calculation. Singular value decomposition can be implemented using numerical computation libraries (such as Python's numpy.linalg.svd function), ensuring the stability and accuracy of the decomposition results. The number of principal components, k, should be determined based on business needs and computational efficiency. For example, while ensuring the cumulative contribution rate, a smaller k value should be chosen to reduce subsequent computational complexity.

[0125] Taking three dimensions (surface displacement, groundwater pressure, and rock formation stress) as an example, the covariance matrix is ​​a 3×3 matrix, and after singular value decomposition, three eigenvalues ​​are obtained: λ1≥λ2≥λ3. If the preset threshold is 85% and the cumulative contribution rate of λ1+λ2 reaches 85%, the first two eigenvectors are selected to form the principal component space, and the three-dimensional fluctuation coefficient vector is reduced to a two-dimensional principal component vector. This dimensionality reduction process can remove redundant information between the original dimensions and retain the main spatiotemporal fluctuation patterns that affect coal mining subsidence. For example, the coordinated change pattern of surface displacement and rock formation stress may constitute the first principal component, while the independent change of groundwater pressure may constitute the second principal component.

[0126] The core function of principal component vector extraction is to transform high-dimensional correlated variables into low-dimensional uncorrelated variables through linear transformation, enabling subsequent spatiotemporal correlation weight calculations to be performed in a simpler, independent feature space. This process, rigorously based on matrix algebra and statistical theory, achieves data dimensionality reduction and feature extraction through spectral decomposition of the covariance matrix, ensuring that the principal component vector effectively characterizes the primary direction of variation of the original spatiotemporal fluctuation coefficients, providing a mathematical foundation for the fusion analysis of multidimensional data.

[0127] Example 4:

[0128] The steps for calculating spatiotemporal association weights are as follows:

[0129] The principal component vector is standardized. Assume that through the principal component vector extraction process of Example 3, a k-dimensional principal component vector pc is obtained (for example, k = 2, corresponding to the comprehensive principal component of surface displacement and rock formation stress and the independent principal component of groundwater pressure). The purpose of the standardization process is to eliminate the dimensional differences between different principal components, so that their mean is 0 and their variance is 1. The specific operation is to calculate the mean and standard deviation of each principal component component, and then subtract the mean from each component value and divide it by the standard deviation. For example, if the two components of the principal component vector are pc1 = 3.2 and pc2 = -1.5, the means of the principal components are 2.8 and -1.2, respectively, and the standard deviations are 0.5 and 0.3, respectively, then the standardized components are (3.2-2.8) / 0.5 = 0.8 and (-1.5+1.2) / 0.3 = -1.0, and the standardized principal component vector [0.8, -1.0] is obtained.

[0130] Calculate the cosine similarity between the normalized principal component vector and the preset reference vector. ref The principal component feature combination representing the typical subsidence pattern is determined by principal component analysis of historical coal mining subsidence data. For example, in historical data, when significant subsidence occurs in the coal mining area, the principal component vector usually presents a pattern of [1.2, -0.8], and this vector is set as the reference vector. Cosine similarity is used to measure the directional similarity of two vectors. It is calculated as the dot product of the two vectors divided by the product of their respective moduli. Assuming that the standardized principal component vector is [0.8, -1.0] and the reference vector is [1.2, -0.8], the dot product is 0.8×1.2+(-1.0)×(-0.8)=0.96+0.8=1.76, and the modulus of the principal component vector is The modulus of the reference vector is The cosine similarity is 1.76 / (1.28×1.44)≈0.95, indicating that the current principal component vector has a high similarity with the typical subsidence pattern.

[0131] The cosine similarity is input into a multilayer perceptron (MLP) for nonlinear mapping to generate initial weights. The MLP consists of an input layer, two hidden layers, and an output layer. The input layer receives a one-dimensional cosine similarity value (e.g., 0.95). The first hidden layer contains 64 neurons and uses the ReLU activation function to map the input value into a nonlinear space. For example, the output value may become max(0, 0.95×W1+b1) (where W1 and b1 are the weights and biases of this layer). The second hidden layer contains 32 neurons and also uses the ReLU activation function. The output layer contains one neuron and uses the hyperbolic tangent (tanh) activation function to compress the output value to the range [-1, 1] to generate the initial weights. For example, if the training parameters of the MLP result in an output of 0.8 for an input of 0.95, the initial weight is 0.8, indicating that the current spatiotemporal fluctuation pattern has a high influence on subsidence prediction.

[0132] The initial weights are smoothed by the gated recurrent unit (GRU). GRU is a recurrent neural network unit that can capture the dependencies in the time series and eliminate the random fluctuations of the initial weights. Assume that the current time step is t and the input is the initial weight w at time t. t =0.8, the hidden state of the previous moment is h t-1 = 0.7. GRU updates the gate z t and reset gate r t Calculate the current hidden state: The update gate determines the degree to which the previous hidden state retains the current state, and the reset gate determines the degree to which the previous hidden state forgets the current candidate state. In the specific calculation, the output of the update gate and the reset gate is the result of the sigmoid function (value range [0, 1]), and the candidate state is generated by the tanh function. For example, if the update gate calculation result is 0.9 and the reset gate calculation result is 0.3, then the previous hidden state h t-1 = 0.7 has a higher degree of preservation, and the candidate state is mainly based on the current input w t = 0.8, the final current hidden state h t = 0.9 × 0.7 + (1 - 0.9) × tanh (...) ≈ 0.63 + 0.1 × 0.8 = 0.71, which is used as the spatiotemporal correlation weight output at time t. GRU's temporal smoothing prevents drastic changes in weight values ​​between adjacent time steps. For example, if the weight was 0.7 at the previous moment, the weight at the current moment will be 0.71 after smoothing. This prevents weight fluctuations caused by noise and enhances the temporal continuity and physical significance of the weights.

[0133] Take a specific application scenario as an example: During the mining process of a coal mine, the spatiotemporal fluctuation coefficient of the surface displacement and the fluctuation coefficient of the rock layer stress show a strong positive correlation (such as increasing at the same time). Through the principal component analysis of the covariance matrix, the first principal component captures the coordinated changes of these two dimensions. The standardized principal component vector is [1.1, -0.9], and the cosine similarity with the preset reference vector [1.0, -1.0] is 0.98. After input into the multi-layer perceptron, the initial weight 0.95 is generated, and then smoothed by GRU. Considering the weight of 0.93 at the previous moment, the update gate output is 0.85, and the reset gate output is 0.2. The current weight is calculated to be This weight indicates that the coordinated fluctuation of surface displacement and rock layer stress has a higher influence weight on subsidence prediction and needs to be given priority consideration in fusion prediction.

[0134] The entire spatiotemporal correlation weight calculation process achieves weight assignment for multi-dimensional spatiotemporal fluctuation coefficients through standardization to eliminate dimension, cosine similarity to measure pattern matching, multi-layer perceptrons to introduce nonlinear mapping, and GRU to smooth temporal fluctuations. This process closely relies on the results of principal component analysis and combines historical subsidence patterns and temporal characteristics to ensure that the weights reasonably reflect the actual contribution of fluctuations in different dimensions to subsidence prediction, providing a scientific basis for dynamically revising the original predicted values.

[0135] Example 5:

[0136] The specific implementation methods of generating fusion prediction values ​​and constructing a three-dimensional subsidence surface model are as follows:

[0137] Taking a mining area of ​​a coal mine as an example, assuming that the spatiotemporal correlation weights of the three dimensions (surface displacement, groundwater pressure, and rock formation stress) obtained by the method of the present invention are 0.6, 0.3, and 0.1 respectively, and the original predicted value is a time series of surface settlement (unit: mm), and the original predicted value F at a certain moment is 50 mm. The spatiotemporal correlation weight is tensor multiplied with the original predicted value, that is, the weight of each dimension is multiplied by the original predicted value component of the corresponding dimension. Since a single dimension (surface settlement) is used as an example here, the weight-corrected predicted value is 50×0.6=30 mm. At this time, the residual between the weight-corrected predicted value and the original predicted value is 30-50=-20 mm. This residual reflects the difference between the weight-corrected predicted value and the original prediction value, and needs to be dynamically compensated by the Kalman filter.

[0138] The Kalman filter's state equation assumes a random walk residual model, meaning the current residual equals the previous residual plus process noise. In the observation equation, the observation is the residual between the weighted predicted value and the original predicted value, including observation noise. Taking time t as an example, assuming the residual estimate at the previous time (t-1) is -18 mm, the process noise covariance is 10, and the observation noise covariance is 5, the Kalman gain is calculated as the process noise covariance divided by (process noise covariance + observation noise covariance), or 10 / (10+5) = 0.667. The predicted residual at time t is the previous residual estimate plus the process noise (assumed to be 0), or -18 mm. The observed value is -20 mm, so the compensated residual estimate is: predicted value + Kalman gain × (observed value - predicted value) = -18 + 0.667 × (-20 + 18) = -18 - 1.334 = -19.334 mm. The compensated residual is added to the original predicted value to generate a fused predicted value of 50+(-19.334)=30.666 mm. This value comprehensively considers the correction of spatiotemporal correlation weights and dynamic compensation of residuals, and is closer to the actual subsidence trend.

[0139] When constructing a three-dimensional subsidence surface model based on the fused prediction value, the fused prediction value of each spatial grid at a certain moment (such as the surface settlement) is mapped to a three-dimensional coordinate system according to the grid coordinates (x, y). The x and y axes are geographic coordinates (unit: meter), and the z axis is the settlement (unit: millimeter), generating discrete point cloud data. For example, if a grid coordinate is (100, 200) and the fused prediction value is 35 mm, the corresponding three-dimensional point is (100, 200, 35). The radial basis function interpolation method is used to fit the surface of the discrete point cloud, and the Gaussian function is selected as the basis function, which is in the form of Φ(r) = exp(-(r / σ) 2 ), where r is the distance between the interpolation point and the data point, and σ is a smoothing parameter (e.g., 50 meters). For the interpolation point (x, y), the fitted value is the sum of the weights of each data point multiplied by the basis function values. The weights are calculated using the least squares method to minimize the error of the fitted surface at the data point.

[0140] After the fitting is completed, the Gaussian curvature and average curvature of the surface are calculated. Gaussian curvature reflects the degree of curvature of the surface in the two main curvature directions, and the average curvature reflects the average curvature of the surface. Taking a certain grid point as an example, if its Gaussian curvature is 0.002 mm-2 and the average curvature is 0.001 mm-1, it indicates that the surface in this area is slightly convex. The curvature distribution is smoothed by the moving least squares method, and a 5×5 spatial window is selected. The curvature value in each window is locally polynomial fitted to eliminate abnormal fluctuations caused by data noise. For example, the original Gaussian curvature calculated value of a point is 0.005 mm-2, but the curvature values ​​of its surrounding neighborhood are all between 0.001 and 0.002 mm-2. After smoothing, the curvature value of the point is corrected to 0.0015 mm-2 to eliminate noise interference.

[0141] Curvature extreme points and gradient break areas are extracted based on smoothed curvature data. Curvature extreme points are local maxima or minima of Gaussian or mean curvature. For example, within a certain area, a point with a Gaussian curvature of 0.008 mm², significantly higher than the 0.001-0.002 mm² of surrounding points, is considered a curvature extreme point. Gradient break areas are determined by calculating the spatial gradient rate of subsidence. A gradient break area is identified if the subsidence gradient of adjacent grid cells within a certain area suddenly increases (e.g., from 1 mm / m to 5 mm / m). When the number of curvature extreme points exceeds a preset threshold (e.g., 10 per square kilometer) or the proportion of gradient break areas exceeds a critical value (e.g., 5%), the area is identified as a high-risk area for subsidence and a graded warning signal is generated. For example, if 15 curvature extreme points per square kilometer are detected within a mining block and the proportion of gradient break areas reaches 7%, a yellow warning is triggered, indicating the need to adjust mining schedules or strengthen monitoring.

[0142] Furthermore, when establishing a joint optimization model combining geological dynamic data and mining schedule parameters, the objective function is defined as the weighted squared error between the predicted subsidence value and the mining schedule, where the mining schedule parameters include mining rate and mining sequence. Constraints are defined as the rock formation stress safety threshold, ensuring that the rock formation stress does not exceed its strength limit during mining. The model is solved using the alternating direction multiplier method (ADMM), decomposing the objective function into a mining schedule subproblem and a subsidence prediction subproblem. The problem is solved iteratively until convergence. For example, the mining schedule is first fixed, and the model parameters that minimize the subsidence prediction error are solved. Then, the model parameters are fixed, and the optimal mining rate and sequence that satisfy the stress constraints are solved. This iterative process continues until both the error and stress conditions are met, resulting in the optimal mining path. The optimal mining path is dynamically coupled with the fused predicted values ​​to generate a mining plan that minimizes subsidence risk. For example, adjusting the mining rate in a certain area from 2 meters per day to 1.5 meters per day can reduce subsidence risk.

[0143] It should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "includes," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that includes a list of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus.

[0144] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. An intelligent prediction method for the spatiotemporal evolution of coal mining subsidence, characterized by: include: Collect multi-dimensional geological dynamic data of the coal mining area, wherein the geological dynamic data includes surface displacement, groundwater pressure and rock formation stress; The GCLSTM algorithm is used to extract spatiotemporal features from geological dynamic data to obtain spatiotemporal feature sequences, and the spatiotemporal fluctuation coefficient is calculated based on the gradient changes of the spatiotemporal feature sequences. The covariance matrix of the spatiotemporal fluctuation coefficients of different dimensions is spectrally decomposed to extract the principal component vectors, which are then nonlinearly mapped using the hyperbolic tangent function to obtain the spatiotemporal correlation weights of each dimension. The original prediction value of the spatiotemporal feature sequence is dynamically modified according to the spatiotemporal correlation weight to generate a fused prediction value, and the spatiotemporal evolution trend of coal mining subsidence is determined based on the fused prediction value; The fusion prediction value F′ satisfies the relationship: F′=F×tanh(SC×W) Where F represents the original prediction value, SC represents the spatiotemporal fluctuation coefficient, W represents the spatiotemporal correlation weight, and tanh() represents the hyperbolic tangent function.

2. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1 is characterized in that: The spatiotemporal feature extraction includes: Intercept geological dynamic data in a preset time window, construct a space-time tensor and perform normalization processing; A gated convolutional layer is used to capture local spatiotemporal features of the spatiotemporal tensor and output the Hadamard product of the gated weights and the convolutional features; The gated convolution result is input into the bidirectional LSTM unit to extract the forward and backward long-term dependency features respectively, and the feature sequence is weightedly fused through the attention mechanism to generate a spatiotemporal feature sequence.

3. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1 is characterized in that: The calculation of the spatiotemporal fluctuation coefficient includes: Select the spatiotemporal feature sequence at any moment as the benchmark feature, and calculate the gradient modulus of the feature vectors at the next n consecutive moments; The ratio of the mean value of the gradient modulus to the Euclidean distance of the benchmark feature is used as the time dimension fluctuation; Calculate the KL divergence between the benchmark feature and the adjacent spatial grid feature, and take its reciprocal as the spatial dimension fluctuation; The geometric mean of the time dimension fluctuation and the space dimension fluctuation is taken as the space-time fluctuation coefficient.

4. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1, characterized in that: The principal component vector extraction includes: Construct the covariance matrix of spatiotemporal fluctuation coefficients of different dimensions, perform singular value decomposition on it and obtain the eigenvalue diagonal matrix; Select eigenvectors whose eigenvalue cumulative contribution rate exceeds a preset threshold to form the principal component space; Project the covariance matrix into the principal component space to obtain the principal component vector after dimensionality reduction.

5. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1 is characterized in that: The spatiotemporal association weight calculation includes: Normalize the principal component vector and calculate its cosine similarity with the preset reference vector; The cosine similarity is input into the multi-layer perceptron, and the initial weight is generated after mapping through the fully connected layer and activation function; The initial weights are temporally smoothed through the gated recurrent unit to output the spatiotemporal correlation weights of each dimension.

6. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1, characterized in that: The generation of the fusion prediction value includes: Perform tensor multiplication of the spatiotemporal correlation weight and the original prediction value to obtain the weight-corrected prediction value; Calculate the residual between the weighted corrected prediction value and the original prediction value, and dynamically compensate the residual through the Kalman filter; The compensated residual is added to the original prediction value to generate a fused prediction value.

7. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1, characterized in that: Also includes: Based on the fusion prediction values, a three-dimensional subsidence surface model is constructed to extract the curvature extreme points and gradient mutation areas of the surface; When the number of extreme curvature points exceeds the preset threshold or the area proportion of the gradient mutation region is greater than the critical value, it is judged as a high-risk area for subsidence and a graded warning signal is generated.

8. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 7, characterized in that: The three-dimensional subsidence surface model construction includes: The fused prediction value is mapped to a three-dimensional coordinate system according to the spatial grid coordinates to generate discrete point cloud data; The radial basis function interpolation method is used to fit the surface of discrete point clouds, and the Gaussian curvature and mean curvature of the fitting surface are calculated; The curvature distribution is smoothed by moving least squares method to eliminate noise interference.

9. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 1, characterized in that: Also includes: A joint optimization model of geological dynamic data and mining progress parameters was established, and the alternating direction multiplier method was used to solve the optimal mining path; The optimal mining path is dynamically coupled with the fused prediction value to generate a mining plan that minimizes the subsidence risk.

10. The method for intelligent prediction of spatiotemporal evolution of coal mining subsidence according to claim 9, characterized in that: The joint optimization model solution includes: The objective function is defined as the weighted square error sum of the subsidence prediction value and the mining progress, and the constraint condition is the rock formation stress safety threshold; Perform Lagrangian dual transformation on the objective function and decompose it into mining progress subproblem and subsidence prediction subproblem; The two sub-problems are solved alternately and iteratively until convergence, and the optimal mining path that meets the constraints is output.

Images