Coal mine surrounding rock deformation prediction method based on LSTM

Abstract

The application relates to the technical field of coal mine surrounding rock deformation monitoring, and discloses a coal mine surrounding rock deformation prediction method based on LSTM. The method comprises the following steps: performing space-time alignment on ground and unmanned aerial vehicle data; constructing an input matrix that fuses multi-source observation, engineering parameters and derived features; outputting a coarse prediction value through double-layer LSTM; establishing a Kalman filter model to fuse the coarse prediction and ground measurement, and calculating a calibration correction amount to obtain a calibrated prediction subsidence amount. The application improves the accuracy, robustness and working condition self-adaptive capability of coal mine surrounding rock deformation prediction.

Description

Technical Field

[0001] This application relates to the field of coal mine surrounding rock deformation monitoring technology, and in particular to a coal mine surrounding rock deformation prediction method based on LSTM. Background Technology

[0002] Predicting coal mine surrounding rock deformation is a crucial technical means for mine safety management and disaster prevention. Existing LSTM-based methods for predicting coal mine surrounding rock deformation construct recurrent neural network models to learn the temporal evolution of surrounding rock deformation using historical monitoring data, thereby predicting future deformation trends. Traditional methods typically use a single data source as model input, such as displacement monitoring data from ground measuring points or surface deformation data acquired by UAVs. These methods are trained using a standard LSTM network architecture, using observations within a historical time window as input features to output predicted deformation values ​​for future moments. Such methods can achieve a certain level of prediction accuracy under conditions of good data quality and stable operating conditions, providing an intelligent technical path for surrounding rock deformation monitoring.

[0003] However, existing LSTM-based methods for predicting coal mine surrounding rock deformation have significant shortcomings. First, the limited prediction accuracy due to a single data source is a prominent issue. While ground-based monitoring has high accuracy, its observation frequency is low and its spatial coverage is limited. Although UAV remote sensing data has wide coverage and a relatively high acquisition frequency, its accuracy is limited by DEM model errors. Using either data source alone makes it difficult to simultaneously meet the temporal resolution and accuracy requirements of prediction. Second, existing methods lack explicit modeling of systematic errors in multi-source data. They simply concatenate observation data of different accuracies and frequencies and input them into the LSTM, failing to effectively handle the spatiotemporal alignment and accuracy differences between data sources. This results in the model's learned feature representation being mixed with real deformation information and data errors, leading to insufficient reliability of prediction results. Third, fixed-parameter LSTM models exhibit poor generalization ability when facing dynamic changes in working conditions during coal mining. The surface deformation mechanism differs in different mining stages: the deformation rate is slow in the initial mining stage, the deformation is severe in the active mining stage, and the deformation gradually stabilizes in the declining mining stage. However, existing methods use the same network structure and hyperparameters to process data from all stages, failing to adaptively adjust the model configuration to match the characteristics of the current working conditions. Summary of the Invention

[0004] This application provides an LSTM-based method for predicting coal mine surrounding rock deformation, which addresses the technical problems of existing LSTM-based methods for predicting coal mine surrounding rock deformation, such as difficulty in balancing temporal resolution and accuracy due to a single data source, lack of explicit modeling of system errors in multi-source data, poor adaptability of fixed parameter models to working conditions, and insufficient samples in scenarios with missing data. By constructing a multi-source data spatiotemporal alignment and feature fusion mechanism, a Kalman filter accuracy calibration module, a parameter pool dynamic optimization strategy, and a soil data reinforcement training method, the accuracy, robustness, and adaptability of coal mine surrounding rock deformation prediction are improved.

[0005] This application provides a method for predicting the deformation of surrounding rock in coal mines based on LSTM, the method comprising:

[0006] Step S1: Perform spatiotemporal alignment processing on the ground subsidence data and UAV elevation change data to obtain a ground subsidence sequence and a UAV elevation change sequence with a unified time reference. Step S2: Construct an input feature matrix containing multi-source observation features, mining engineering features, and spatiotemporal derived features based on the ground subsidence sequence and the UAV elevation change sequence; Step S3: Input the UAV data stream in the input feature matrix into the coarse-grained LSTM path, perform forward propagation calculation through the forget gate, input gate and output gate of two layers of LSTM units, and output the coarse predicted sinking vector for the future time period through the fully connected layer. Step S4: Establish a Kalman filter state space model, using the coarse predicted subsidence vector as the medium-precision observation channel and the measured ground subsidence value as the high-precision observation channel. Calculate the Kalman gain using a joint observation matrix and a piecewise filtering update strategy, fuse the data from the two channels to obtain the optimal state estimate, calculate the deviation between the optimal state estimate and the coarse predicted subsidence vector as the calibration correction, and superimpose the calibration correction onto the coarse predicted subsidence vector to obtain the calibrated predicted subsidence.

[0007] The technical solution provided in this application obtains two types of sequences with a unified time reference by performing spatiotemporal alignment processing on ground subsidence data and UAV elevation change data. This solves the problem of data time dimension mismatch caused by differences in sampling frequencies between different observation methods in existing technologies, enabling high-precision low-frequency ground observation data and medium-precision high-frequency UAV observation data to be synchronized on the time axis. Based on the aligned ground subsidence sequence and UAV elevation change sequence, an input feature matrix containing multi-source observation features, mining engineering features, and spatiotemporally derived features is constructed. The systematic error between the data sources is explicitly expressed by calculating the difference features between the two types of observation data. Subsidence rate features and crack development growth rate features are extracted. By incorporating the dynamic evolution information of the deformation process, the LSTM model can not only learn the temporal patterns of a single data source, but also capture the complementary relationships and error coupling patterns between multiple data sources. This overcomes the problem of chaotic feature representation caused by simply splicing multiple data sources in existing technologies. The UAV data stream in the input feature matrix is ​​input into a coarse-grained LSTM path and forward propagated through the forget gate, input gate, and output gate of two layers of LSTM units. By using the gating mechanism to selectively memorize long-term dependencies, the model can extract the macroscopic trend and periodic features of subsidence evolution from high-frequency observations of UAV data. The coarse predicted subsidence vector for future periods is output through a fully connected layer, providing a preliminary prediction benchmark for subsequent accuracy calibration.

[0008] The core innovation of this application lies in establishing a Kalman filter state-space model. The coarse predicted subsidence vector is used as the medium-precision observation channel, and the measured ground subsidence value is used as the high-precision observation channel. The optimal state estimate is obtained by fusing the data from the two channels through a joint observation matrix and a piecewise filtering update strategy. This design fully leverages the theoretical advantages of Kalman filtering in fusing observation data of different accuracies. It automatically allocates fusion weights based on the differences in the observation noise covariance of each channel, ensuring that high-precision ground observation data dominates state updates when available, and that medium-precision coarse predicted data maintains state tracking when ground observations are missing. This achieves dynamic complementarity between different data sources in the time dimension. The deviation between the optimal state estimate and the coarse predicted subsidence vector is calculated as a calibration correction and superimposed on the coarse predicted vector to obtain the calibrated predicted subsidence. The calibration mechanism is not a simple post-processing correction, but rather establishes an explicit transmission relationship between prediction and observation errors based on a state-space model. This allows the calibration correction to reflect systematic biases and random disturbances between data sources, overcoming the accuracy limitations caused by directly using LSTM output as the prediction result in existing technologies. Through a progressive technology chain of multi-source data spatiotemporal alignment, difference feature construction, dual-path LSTM processing, and Kalman filter calibration, this application fully leverages the temporal modeling capability of the LSTM algorithm and the data fusion capability of the Kalman filter algorithm in the specific application field of coal mine surrounding rock deformation prediction. This enables the prediction results to simultaneously possess the high temporal resolution of UAV data and the high precision of ground observation data, providing an effective technical solution for refined monitoring of surface subsidence under conditions of insufficient mining. Attached Figure Description

[0009] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0010] Figure 1 This is a schematic diagram of an embodiment of the LSTM-based method for predicting deformation of surrounding rock in coal mines in this application. Figure 2 This is a schematic diagram comparing the prediction accuracy before and after parameter optimization in an embodiment of this application; Figure 3 This is a schematic diagram comparing the prediction accuracy at different mining stages in the embodiments of this application. Detailed Implementation

[0011] This application provides an LSTM-based method for predicting deformation of surrounding rock in coal mines. The terms "first," "second," "third," "fourth," etc. (if present) in the specification, claims, and accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in a sequence other than that illustrated or described herein. Furthermore, the terms "comprising" or "having" and any variations thereof are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0012] For ease of understanding, the specific process of the embodiments of this application is described below. Please refer to [link / reference]. Figure 1 One embodiment of the LSTM-based method for predicting deformation of surrounding rock in coal mines in this application includes: Step S1: Perform spatiotemporal alignment processing on the ground subsidence data and UAV elevation change data to obtain a ground subsidence sequence and a UAV elevation change sequence with a unified time reference. Among them, the problem of inconsistent sampling frequencies of multi-source data is solved by spatiotemporal alignment processing. High-precision subsidence data is obtained by observing ground measuring points every 7 days, and DEM data is obtained by UAVs every 10 days. The reference time step is set to 1 day. The ground measuring point data is extended into a daily sequence by cubic spline interpolation, and the UAV DEM data is extended into a daily grid by bilinear interpolation. At the same time, the DEM grid nodes within a 5-meter radius around the ground measuring points are extracted, and the elevation change of each node is weighted by the inverse square of the distance as the weight. This achieves spatial alignment of the area DEM to the point measuring points, so that data of different precision and frequency are consistent in both time and space dimensions.

[0013] Step S2: Construct an input feature matrix containing multi-source observation features, mining engineering features, and spatiotemporal derived features based on the ground subsidence sequence and the UAV elevation change sequence; Specifically, an input feature matrix integrating multi-scale information is constructed. The difference between ground subsidence and UAV elevation change is calculated as a multi-source data difference feature to explicitly express the systematic error. The ratio of subsidence change to time step length in adjacent time steps is calculated as a subsidence rate feature to reflect deformation dynamics. The ratio of crack density change to time interval in 10-day intervals is calculated as a crack development growth rate to quantify the degree of surface damage. Combined with engineering parameters such as working face advance speed, coal seam thickness, working face width, and mining depth, the matrix is ​​organized into a matrix with 13 feature dimensions according to a 30-day time window to provide structured input for LSTM.

[0014] Step S3: Input the UAV data stream in the input feature matrix into the coarse-grained LSTM path, perform forward propagation calculation through the forget gate, input gate and output gate of two layers of LSTM units, and output the coarse predicted sinking vector for the future time period through the fully connected layer. Specifically, the elevation change of the UAV, crack density, and crack development growth rate in the input feature matrix are extracted as coarse-grained input submatrices. The input consists of two-layer LSTM networks with 64 and 32 hidden units, respectively. At each time step, the first-layer LSTM unit controls the retention of the previous time step through a forget gate, controls the reception of the current input through an input gate, updates the unit state, and controls the output of the hidden state through an output gate. The second-layer LSTM repeats this process with the output of the first layer as input. Finally, the hidden state at the last time step is passed through a fully connected layer and a ReLU activation function to output a coarse prediction vector of the subsidence amount for the next 5 days. This coarse prediction inherits the advantage of high frequency of UAV data, but its accuracy is limited by a DEM error of 5 cm.

[0015] Step S4: Establish a Kalman filter state space model, using the coarse predicted subsidence vector as the medium-precision observation channel and the actual ground subsidence value as the high-precision observation channel. Calculate the Kalman gain through a joint observation matrix and a piecewise filtering update strategy, fuse the data from the two channels to obtain the optimal state estimate, calculate the deviation between the optimal state estimate and the coarse predicted subsidence vector as the calibration correction, and superimpose the calibration correction onto the coarse predicted subsidence vector to obtain the calibrated predicted subsidence.

[0016] Specifically, a Kalman filter state-space model is established to fuse observation data of different accuracies. On days 1 to 6 when no ground-measured data is available, only coarse prediction values ​​are used to update the state estimate through single-channel Kalman gain. On day 7 when ground-measured data is available, a joint observation matrix including the ground observation matrix and the coarse prediction observation matrix is ​​constructed. Based on the ground observation noise covariance and the coarse prediction observation noise covariance, a diagonal block joint noise covariance matrix is ​​constructed. The optimal state estimate is obtained by fusing the prior state estimate, the ground-measured subsidence value, and the coarse prediction subsidence amount using dual-channel Kalman gain. The deviation between the optimal estimate and the coarse prediction value is the calibration correction amount. The correction amount is superimposed on the coarse prediction vector to obtain the calibrated predicted subsidence amount, which reduces the prediction accuracy from 4.8 cm to 2.1 cm.

[0017] In one specific embodiment, step S1 includes: The ground subsidence data at the measuring points were extended by cubic spline interpolation according to the reference time step to obtain the daily ground subsidence data sequence. The UAV elevation change data is extended by bilinear interpolation according to the reference time step to obtain the daily UAV elevation change data sequence. Extract DEM grid nodes within a set radius around the plane coordinates of the ground measuring point, calculate the Euclidean distance from each grid node to the measuring point, and use the inverse square of the distance as the weight to perform a weighted average of the elevation changes of the grid nodes to obtain the UAV elevation change sequence corresponding to the measuring point location. Align the daily ground subsidence data series and the UAV elevation change series in the time dimension to obtain the ground subsidence series and the UAV elevation change series.

[0018] Specifically, cubic spline interpolation extension achieves smooth data transition by constructing a cubic polynomial function between known ground measurement point observation times. A cubic polynomial segment is formed between every two adjacent observation points, and the function value, first derivative, and second derivative are kept continuous at the connection points of each segment. This ensures that the interpolated daily ground subsidence data sequence is not only numerically accurate but also has a smooth trend, avoiding the broken-line abrupt changes caused by linear interpolation. Bilinear interpolation extension is performed on UAV DEM raster data in the time dimension. For the time when interpolation is required, DEM data obtained from two observations before and after that time are selected, and the elevation change at each grid location is weighted according to the time distance to obtain the interpolation result, so that the daily UAV elevation change data sequence presents a continuous distribution on the time axis.

[0019] Spatial mapping converts area DEM data into point data corresponding to the measurement point location using a distance-inverse square weighted method. For a given ground measurement point, there are multiple DEM grid nodes within a defined radius. The influence of the elevation change of each grid node on the measurement point is inversely proportional to the square of the distance; the closer the grid node, the greater its weight. The UAV elevation change at the measurement point is calculated by weighted averaging of the elevation changes of all grid nodes. For example, if there are 12 DEM grid nodes within a 5-meter radius around a measurement point, with distances ranging from 1 meter to 4.8 meters, after weighting according to the distance-inverse square weighted method, the node at 1 meter has a weight of 1.0, while the node at 4.8 meters has a weight of only 0.043. After weighted averaging, a UAV elevation change sequence precisely corresponding to the measurement point location is obtained. This is combined with the daily ground subsidence data sequence to achieve time dimension alignment, forming two types of observation data with a unified time reference.

[0020] In one specific embodiment, step S2 includes: The difference between the ground subsidence sequence and the UAV elevation change sequence is calculated to obtain the multi-source data difference feature sequence; The ratio of the change in ground subsidence between adjacent time steps to the time step length is calculated to obtain the subsidence rate characteristic sequence. The ratio of the change in crack density over a set time interval to the time interval is calculated to obtain the crack development growth rate characteristics. The ground subsidence sequence, UAV elevation change sequence, multi-source data difference feature sequence, subsidence rate feature sequence, working face advance speed, coal seam mining thickness, working face width, mining depth and fracture development growth rate features are organized into a matrix form according to a set time window length to obtain the input feature matrix.

[0021] Specifically, the multi-source data difference feature sequence is obtained by calculating the difference between ground subsidence and UAV elevation change at each time step. This difference reflects the systematic bias of the two observation methods. Ground measurement points use high-precision electronic levels to observe subsidence with millimeter-level accuracy, while the elevation error of the UAV DEM model is about 5 centimeters. The difference includes both measurement errors from the UAV data and errors from spatial interpolation. This difference feature, as an explicit input, helps the LSTM model learn and compensate for the systematic bias patterns between multi-source data. The subsidence rate feature sequence is obtained by calculating the change in ground subsidence between adjacent time steps. Dividing by the time step, since the time step is uniformly 1 day, the subsidence rate is the daily subsidence increment. This feature reflects the dynamic evolution speed of surface deformation. The subsidence rate is larger during the rapid advancement of the working face and gradually decreases during the decline of mining impact. The crack development growth rate feature is obtained by calculating the change in crack density over a 10-day time interval and dividing it by the time interval. Crack density is extracted from high-definition UAV images using the U-Net++ deep learning model, and the unit is cracks per square meter. The crack development growth rate quantifies the development speed of surface damage and has a strong correlation with working face advancement and subsidence deformation.

[0022] The input feature matrix is ​​organized using a sliding 30-day time window. For the prediction time t, various feature data from the past 30 days (t-29 to t) are extracted. The ground subsidence sequence, UAV elevation change sequence, multi-source data difference feature sequence, and subsidence rate feature sequence are arranged according to the time dimension to form a 30-row time series data. At the same time, the working face advance speed, coal seam mining thickness, working face width, mining depth, and fracture development growth rate features are filled in at the corresponding time points, finally forming a 30-row, 13-column input feature matrix, where 30 represents the number of time steps and 13 represents the number of feature dimensions. For example, if a working face needs to predict the subsidence amount for the next 5 days on day 100, then historical data from day 71 to day 100 (30 days in total) is extracted. Each day includes 13 feature values ​​such as ground subsidence amount of 0.012 meters, UAV elevation change of 0.015 meters, difference of 0.003 meters, subsidence rate of 0.0008 meters per day, crack density of 0.6 cracks per square meter, crack growth rate of 0.05 cracks per square meter per day, working face advance speed of 2.5 meters per day, coal seam thickness of 3.2 meters, working face width of 240 meters, and mining depth of 450 meters. These features are organized into a 30x13 matrix and input into the LSTM network.

[0023] In one specific embodiment, the coarse-grained LSTM path includes a first layer LSTM unit and a second layer LSTM unit. The first layer LSTM unit contains a first predetermined number of hidden units, and the second layer LSTM unit contains a second predetermined number of hidden units. The first predetermined number is greater than the second predetermined number.

[0024] Specifically, the coarse-grained LSTM path employs a two-layer stacked structure to achieve deep feature extraction from the UAV data stream. The first LSTM unit contains 64 hidden units, and the second LSTM unit contains 32 hidden units, presenting a progressively decreasing network architecture. The first LSTM unit is responsible for extracting primary temporal features from the input UAV elevation change, crack density, and crack development growth rate. The 64 hidden units provide a large feature representation space, capable of capturing diverse patterns in the subsidence evolution process, including periodic fluctuations, trend changes, and abrupt changes. The second LSTM unit receives the 64-dimensional hidden state output from the first layer as input, and extracts further features through the 32 hidden units. The hidden unit further refines and compresses the feature representation, removes redundant information and focuses on the most critical high-order abstract features for the prediction task. This layer-by-layer decreasing design follows the feature pyramid principle, with the bottom layer retaining more detailed information and the top layer extracting more compact semantic information. For example, the first layer LSTM may learn a local pattern in which the crack density fluctuates slightly every 3 days. The second layer LSTM integrates these local patterns into a global correlation between the working face advancement cycle and surface damage. Finally, the 32-dimensional hidden state output by the second layer LSTM at the 30th time step contains the core features of the subsidence evolution over the past 30 days. After being mapped by a fully connected layer, it becomes a coarse prediction vector of subsidence for the next 5 days.

[0025] In one specific embodiment, forward propagation computation is performed using the forget gate, input gate, and output gate of two LSTM units, including: The UAV elevation change, crack density, and crack development growth rate in the input feature matrix are extracted as coarse-grained input submatrices. The coarse-grained input submatrix is ​​progressively input into the first-layer LSTM unit along the time dimension. At each time step, the forget gate activation value, input gate activation value, candidate memory unit value, and output gate activation value are calculated. The current time step unit state is updated based on the sum of the element-wise product of the forget gate activation value and the unit state of the previous time step, and the element-wise product of the input gate activation value and the candidate memory unit value. The current time step hidden state is calculated based on the element-wise product of the output gate activation value and the hyperbolic tangent of the current time step unit state. The hidden states output by the first-layer LSTM unit at each time step are used as the input of the second-layer LSTM unit. Following the calculation process of the first-layer LSTM unit, forward propagation is performed in the second-layer LSTM unit to obtain the output hidden states of the second-layer LSTM unit at the last time step. The output hidden state of the second LSTM unit is input into the fully connected layer, and a nonlinear transformation is performed through the ReLU activation function to obtain the coarse prediction sink vector.

[0026] Specifically, the extraction of the coarse-grained input submatrix targets the feature dimensions related to UAV observations in the input feature matrix. From the complete 30-row, 13-column input feature matrix, three feature columns are extracted: UAV elevation change, crack density, and crack development growth rate, forming a 30-row, 5-column coarse-grained input submatrix. The five columns include one column for UAV elevation change, one column for crack density, one column for crack development growth rate, and two columns for the spatial correlation features of these three features with adjacent measuring points. This submatrix is ​​specifically used for coarse-grained LSTM path processing, reflecting the design idea of ​​modeling different data source characteristics separately in the dual-path network. Although UAV data has limited accuracy, it has the advantages of high frequency and wide coverage, making it suitable for capturing the macroscopic trend and spatial distribution characteristics of subsidence evolution.

[0027] The forward propagation of the first-layer LSTM unit processes each row of the coarse-grained input submatrix sequentially according to the time step order. At the first time step, the first row of the 5-dimensional feature vector is input. The forget gate uses a sigmoid activation function to calculate an activation value between 0 and 1, controlling the retention ratio of the unit state from the previous time step. The input gate also uses a sigmoid activation function to calculate an activation value, controlling the reception ratio of the current input information. The candidate memory unit uses a hyperbolic tangent activation function to calculate candidate values ​​between -1 and +1, representing new information at the current time step. The output gate uses a sigmoid activation function to calculate the output ratio from the control state to the hidden state. The unit state update at the current time step retains some historical memory by element-wise multiplying the forget gate activation value by the unit state from the previous time step. The input gate activation value is multiplied element-wise with the candidate memory cell value and added to the current new information. The sum of the two parts yields the updated cell state. The hidden state at the current time step is obtained by multiplying the output gate activation value element-wise with the hyperbolic tangent of the current cell state. This hidden state is a 64-dimensional vector containing the accumulated information from the first time step to the current time step. For example, at the 10th time step, the forget gate activation value of 0.85 indicates that 85% of the cell state from the previous time step is retained, the input gate activation value of 0.65 indicates that 65% of the current input is received, and the output gate activation value of 0.78 indicates that 78% of the information of the current cell state is output. This gating mechanism enables LSTM to selectively memorize long-term dependencies. After processing 30 time steps, the first layer of LSTM cells outputs 30 hidden state vectors.

[0028] The second LSTM unit receives 30 hidden state vectors from the first LSTM unit as its input sequence. Each hidden state is a 64-dimensional vector. The second LSTM unit also performs forward propagation sequentially according to time steps, calculating the forget gate activation value, input gate activation value, candidate memory unit value, and output gate activation value at each time step. It updates the unit state and calculates the hidden state using the same gating mechanism. The difference is that the second LSTM unit contains 32 hidden units, so the output hidden state is a 32-dimensional vector. The input to the second layer is the high-level feature representation extracted by the first layer, rather than the original observation data, enabling more abstract feature combinations and pattern recognition. After processing 30 time steps, the output hidden state of the 30th and final time step is taken as the feature representation of the entire coarse-grained LSTM path. This 32-dimensional vector condenses the past 30 days. The core information most relevant to future subsidence prediction from UAV observation data is input into a fully connected layer containing 16 neurons. After linear transformation using weight matrices and bias terms, the ReLU activation function is used to set negative values ​​to zero and retain positive values, achieving nonlinear mapping. The 16-dimensional vector output by the fully connected layer is then mapped to a 5-dimensional vector through an output layer. Each dimension corresponds to the coarse predicted subsidence amount from day 1 to day 5. For example, after processing by LSTM, the coarse predicted subsidence vector output for a certain measuring point is 0.018 meters on day 1, 0.022 meters on day 2, 0.025 meters on day 3, 0.027 meters on day 4, and 0.029 meters on day 5, reflecting the increasing trend of subsidence over time. This coarse prediction result inherits the temporal resolution advantage of UAV data, but its accuracy is limited by DEM model errors and requires further correction in subsequent calibration steps.

[0029] Figure 2 This is a schematic diagram comparing the prediction accuracy before and after parameter optimization in the embodiments of this application. Figure 2 The comparison between the LSTM model prediction curve and the measured ground values ​​before and after parameter optimization is shown. The horizontal axis represents the prediction time step for the next 5 days, and the vertical axis represents the subsidence. The measured ground values ​​show an increasing trend from 0.052 meters to 0.074 meters. Before optimization, there was a systematic deviation between the prediction curve and the measured values, with the prediction error reaching 0.3 cm on the 3rd day. After optimization, the prediction curve, through Kalman filter calibration and parameter adaptive adjustment, reduced the prediction error to 0.1 cm on the 3rd day. The overall prediction trajectory highly matches the measured curve, verifying the effectiveness of the multi-source data fusion calibration mechanism and parameter dynamic optimization strategy in this invention. This reduces the prediction accuracy of the LSTM model under non-fully mined conditions from a coarse prediction of 4.8 cm to a final 1.5 cm, meeting the engineering requirements for refined monitoring of coal mine surrounding rock deformation.

[0030] In one specific embodiment, the Kalman filter state-space model includes a state equation and an observation equation. The state equation describes the evolution of the actual subsidence over time, and the observation equation includes a high-precision ground observation channel and a medium-precision coarse prediction observation channel.

[0031] Specifically, the Kalman filter state-space model describes the data fusion process through state equations and observation equations. The state equations represent the evolution of the actual subsidence between adjacent time points; the actual subsidence at the current time point equals the actual subsidence at the previous time point plus process noise. The process noise follows a normal distribution with a mean of 0, and the covariance is 0.0005 square meters to characterize the random disturbance in the surface subsidence process. The observation equations include two channels: a high-precision ground observation channel correlates the electronic level measurement with the actual subsidence, and the observation noise covariance is set to 5 x 10^-7 square meters based on the instrument accuracy. (The last sentence appears to be incomplete and possibly refers to a coarse prediction.) The precision observation channel correlates the coarse prediction value output by LSTM with the actual subsidence. The observation noise covariance is set to 0.0023 square meters based on the coarse prediction error. The ratio of the noise covariance of the two channels is about 4600 times, reflecting that the ground observation precision is much higher than the coarse prediction precision. Kalman filtering automatically allocates fusion weights based on this precision difference. When both ground observation and coarse prediction data are available on day 7, the ground observation channel receives about 98% of the weight. When only coarse prediction data is available on days 1 to 6, the state estimation is updated entirely by relying on the coarse prediction channel. This segmented fusion strategy achieves the complementary advantages of different precision data sources.

[0032] In one specific embodiment, the Kalman gain is calculated using a joint observation matrix and a piecewise filtering update strategy, including: During the period without ground-based measured data, the prior error covariance is calculated based on the state transition matrix and the posterior error covariance of the previous time step. The single-channel Kalman gain is calculated based on the prior error covariance, the coarse prediction observation matrix, and the coarse prediction observation noise covariance. The prior state estimate and the coarse prediction sinking vector are fused based on the single-channel Kalman gain to obtain the updated posterior state estimate and posterior error covariance for the single channel. Within the time period where ground-measured data is available, a joint observation matrix containing the ground observation matrix and the coarse prediction observation matrix is ​​constructed. A diagonal block joint observation noise covariance matrix containing the ground observation noise covariance and the coarse prediction observation noise covariance is also constructed. The dual-channel Kalman gain is calculated based on the prior error covariance, the joint observation matrix, and the joint observation noise covariance matrix. The prior state estimate, the ground-measured subsidence value, and the coarse prediction subsidence vector are fused based on the dual-channel Kalman gain to obtain the updated dual-channel posterior state estimate and posterior error covariance.

[0033] Specifically, during the period from day 1 to day 6 when no ground-based measured data was available, the Kalman filter performed a single-channel update process. First, the prior error covariance was calculated based on the state transition matrix and the posterior error covariance of the previous time step. Since the state transition matrix was set to 1, the prior error covariance equals the posterior error covariance of the previous time step plus the process noise covariance of 0.0005 square meters. This prior error covariance reflects the uncertainty of the prediction based on the state equation. Then, the single-channel Kalman gain was calculated based on the prior error covariance, the coarse prediction observation matrix, and the coarse prediction observation noise covariance of 0.0023 square meters. The Kalman gain is a weighting coefficient between 0 and 1, representing the relative weighting of the coarse prediction observation during state updates. Regarding the reliability of the prior state estimate, when the prior error covariance is large, a Kalman gain close to 1 indicates greater confidence in the observed value; when the observation noise covariance is large, a Kalman gain close to 0 indicates greater confidence in the prior estimate. After calculating the Kalman gain, the prior state estimate and the coarse predicted sinking vector are weighted and fused. The fusion formula is that the posterior state estimate equals the prior state estimate plus the Kalman gain multiplied by the difference between the coarse predicted value and the prior estimate. At the same time, the posterior error covariance is updated to reflect the uncertainty of the fused estimate. For example, on the 3rd day, the prior error covariance is 0.0028 square meters, and the calculated Kalman gain is approximately 0.55, indicating that the coarse predicted observation contributes 55% to the posterior estimate while the prior estimate contributes 45%.

[0034] During the 7th day period when ground-measured data is available, the Kalman filter performs a dual-channel joint update process. A 2x1 matrix is ​​constructed, containing the ground observation matrix and the coarse prediction observation matrix. Two rows of elements are both 1, indicating that both observation channels directly observe the actual subsidence. A 2x2 diagonal matrix is ​​also constructed for the diagonal block joint observation noise covariance matrix. The first element on the diagonal is the ground observation noise covariance (5 x 10^-7 square meters), the second element is the coarse prediction observation noise covariance (0.0023 square meters), and off-diagonal elements are 0, indicating that the noise from the two observation channels is independent. The dual-channel Kalman gain is calculated based on the prior error covariance, the joint observation matrix, and the joint observation noise covariance matrix. Since the accuracy of ground observation is much higher than that of coarse prediction, the calculated Kalman gain vector corresponds to the ground observation data. The magnitude is approximately 0.98, while the corresponding coarse prediction component is approximately 0.02. During fusion, the posterior state estimate equals the prior state estimate plus the inner product of the Kalman gain vector and the difference vector between the observation vector and the prior estimate. In practice, the measured ground subsidence value dominates the state update, while the coarse prediction value only plays a minor auxiliary correction role. For example, on day 7, the measured ground value is 0.052 meters, the coarse prediction value is 0.048 meters, and the prior estimate is 0.050 meters. After fusion, the posterior state estimate is 0.052 meters, almost entirely based on the measured ground value. After the posterior error covariance is updated, it is significantly reduced to a level close to the ground observation noise covariance, reflecting that high-precision observation significantly reduces estimation uncertainty. This dual-channel fusion mechanism ensures that when high-precision observation is available, its accuracy advantage is fully utilized, while when high-precision observation is missing, state tracking can still be maintained by relying on coarse prediction.

[0035] Figure 3 This is a schematic diagram comparing the prediction accuracy at different mining stages in the embodiments of this application. Figure 3 This paper compares the prediction errors of fixed-parameter LSTM and parameter-adaptive LSTM at different mining stages. The horizontal axis represents the three stages: initial mining, active mining, and declining mining, while the vertical axis represents the prediction error RMSE. The prediction errors of fixed-parameter LSTM in the three stages are 3.2 cm, 5.8 cm, and 4.1 cm, respectively. After dynamically optimizing the hyperparameters using the PSO algorithm, the prediction errors of parameter-adaptive LSTM are reduced to 2.1 cm, 2.8 cm, and 2.3 cm, respectively, representing reductions of 34%, 52%, and 44%. This verifies the adaptability of the parameter pool dynamic selection mechanism in this invention to different working conditions. In particular, under the conditions of rapid face advance and severe surface deformation during the active mining stage, the parameter-adaptive strategy significantly improves the model's prediction accuracy by adjusting hyperparameters such as the number of hidden layer units and the learning rate of the LSTM in real time.

[0036] In one specific embodiment, the optimal state estimate is obtained, and the deviation between the optimal state estimate and the coarsely predicted subsidence vector is calculated as a calibration correction amount, including: The posterior state estimate after single-channel update and the posterior state estimate after dual-channel update are used as the optimal state estimate; The difference between the optimal state estimate and the coarse prediction of the subsidence vector at the corresponding time is calculated to obtain the calibration correction amount; Record the calibration correction values ​​of each measuring point within a continuously set time window, construct a spatiotemporal error transfer matrix, perform principal component analysis on the spatiotemporal error transfer matrix, extract the principal components whose cumulative variance contribution rate reaches a set threshold, and obtain the spatial correlation pattern of error transfer.

[0037] Specifically, the posterior state estimate after single-channel update corresponds to the optimal estimate obtained by fusing prior estimates and coarse predictions through Kalman filtering during the time when there is no ground measured data from day 1 to day 6. The posterior state estimate after dual-channel update corresponds to the optimal estimate obtained by jointly fusing prior estimates, ground observations, and coarse predictions during the time when ground measured data is available on day 7. Together, they constitute the optimal state estimate sequence within the complete observation period. This sequence represents the best estimate of the actual subsidence after comprehensively utilizing all available information. Compared with using coarse predictions or ground observations alone, the optimal state estimate fully integrates the advantages of different data sources, maintains continuity during periods when ground observations are missing, and achieves high-precision calibration when ground observations are available.

[0038] The calibration correction is obtained by calculating the difference between the optimal state estimate and the coarse predicted subsidence vector at the corresponding time. This difference quantifies the systematic deviation of the coarse prediction relative to the fused optimal estimate. For example, if the optimal state estimate on day 3 is 0.031 meters while the coarse prediction is 0.028 meters, the calibration correction is positive 0.003 meters, indicating that the coarse prediction underestimates the actual subsidence. If the optimal state estimate on day 7 is 0.052 meters while the coarse prediction is 0.048 meters, the calibration correction is positive 0.004 meters. Adding this correction to the coarse prediction value for future periods can achieve accuracy compensation. Assuming the coarse prediction value on day 8 is 0.055 meters, after adding the correction of 0.004 meters from day 7, the calibrated prediction value is 0.059 meters, which is closer to the actual subsidence trend.

[0039] The spatiotemporal error propagation matrix is ​​constructed by recording the calibration correction amounts at each measuring point within a continuous 30-day time window. This matrix has 100 rows (equal to the number of measuring points) and 30 columns (equal to the number of time steps). The element in the i-th row and j-th column represents the calibration correction amount at the i-th measuring point on day j. This matrix contains both the evolution of correction amounts in the time dimension and the distribution differences of correction amounts at different measuring points in the spatial dimension. Principal component analysis (PCA) is performed on this matrix to extract the main change patterns. PCA projects the 100-dimensional measuring point space onto a low-dimensional principal component space. Each principal component represents a spatial correlation pattern of error propagation. For example, the first principal component might reflect the distribution characteristic of generally larger correction amounts at measuring points directly below the working surface and smaller correction amounts at edge measuring points; the second principal component might reflect the direction of movement. By analyzing the differential propagation patterns of the correction amount in the tendency direction, the top three principal components with a cumulative variance contribution rate of 85% were extracted. These three principal components condense 85% of the information in the spatiotemporal error propagation matrix. The principal component loading vector describes the projection weight of each measuring point on the principal component. High-weight measuring points correspond to the main influence areas of the error pattern. By analyzing the principal component loading vector, the spatial propagation path and influence range of the coarse prediction system error can be identified, providing prior knowledge of error characteristics for subsequent parameter optimization and model improvement. For example, if the first principal component shows that the correction amount of the measuring points within 50 meters in front of the working face advance direction is systematically large, it indicates that the coarse prediction model does not adequately predict the impact of mining on this area, and it is necessary to adjust the model parameters or increase the feature weight of this area.

[0040] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for predicting deformation of surrounding rock in coal mines based on LSTM, characterized in that, The method includes: Step S1: Perform spatiotemporal alignment processing on the ground subsidence data and UAV elevation change data to obtain a ground subsidence sequence and a UAV elevation change sequence with a unified time reference. Step S2: Construct an input feature matrix containing multi-source observation features, mining engineering features, and spatiotemporal derived features based on the ground subsidence sequence and the UAV elevation change sequence; Step S3: Input the UAV data stream in the input feature matrix into the coarse-grained LSTM path, perform forward propagation calculation through the forget gate, input gate and output gate of two layers of LSTM units, and output the coarse predicted sinking vector for the future time period through the fully connected layer. Step S4: Establish a Kalman filter state space model, using the coarse predicted subsidence vector as the medium-precision observation channel and the measured ground subsidence value as the high-precision observation channel. Calculate the Kalman gain using a joint observation matrix and a piecewise filtering update strategy, fuse the data from the two channels to obtain the optimal state estimate, calculate the deviation between the optimal state estimate and the coarse predicted subsidence vector as the calibration correction, and superimpose the calibration correction onto the coarse predicted subsidence vector to obtain the calibrated predicted subsidence.

2. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 1, characterized in that, Step S1 includes: The ground subsidence data at the measuring points were extended by cubic spline interpolation according to the reference time step to obtain the daily ground subsidence data sequence. The UAV elevation change data is extended by bilinear interpolation according to the aforementioned reference time step to obtain a daily UAV elevation change data sequence. Extract DEM grid nodes within a set radius around the plane coordinates of the ground measuring point, calculate the Euclidean distance from each grid node to the measuring point, and use the inverse square of the distance as the weight to perform a weighted average of the elevation changes of the grid nodes to obtain the UAV elevation change sequence corresponding to the measuring point location. Align the daily ground subsidence data sequence and the UAV elevation change sequence in the time dimension to obtain the ground subsidence sequence and the UAV elevation change sequence.

3. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 1, characterized in that, Step S2 includes: The difference between the ground subsidence sequence and the UAV elevation change sequence is calculated to obtain a multi-source data difference feature sequence; The ratio of the change in ground subsidence sequence between adjacent time steps to the time step length is calculated to obtain the subsidence rate characteristic sequence. The crack development growth rate characteristic is obtained by calculating the ratio of the change in crack density over a set time interval to the time interval. The ground subsidence sequence, the UAV elevation change sequence, the multi-source data difference feature sequence, the subsidence rate feature sequence, the working face advance speed, coal seam mining thickness, working face width, mining depth, and the fracture development growth rate feature are organized into a matrix according to a set time window length to obtain the input feature matrix.

4. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 1, characterized in that, In step S3, the coarse-grained LSTM path includes a first layer LSTM unit and a second layer LSTM unit. The first layer LSTM unit contains a first set number of hidden units, and the second layer LSTM unit contains a second set number of hidden units. The first set number is greater than the second set number.

5. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 4, characterized in that, The forward propagation calculation using the forget gate, input gate, and output gate of two-layer LSTM units includes: The UAV elevation change, crack density, and crack development growth rate in the input feature matrix are extracted into a coarse-grained input submatrix. The coarse-grained input submatrix is ​​progressively input into the first layer LSTM unit along the time dimension. At each time step, the forget gate activation value, input gate activation value, candidate memory unit value, and output gate activation value are calculated. The current time step unit state is updated based on the sum of the element-wise product of the forget gate activation value and the unit state of the previous time step, and the element-wise product of the input gate activation value and the candidate memory unit value. The current time step hidden state is calculated based on the element-wise product of the output gate activation value and the hyperbolic tangent of the current time step unit state. The hidden state output by the first LSTM unit at each time step is used as the input of the second LSTM unit. The forward propagation is performed in the second LSTM unit according to the calculation process of the first LSTM unit to obtain the output hidden state of the second LSTM unit at the last time step. The output hidden state of the second LSTM unit is input into the fully connected layer, and a nonlinear transformation is performed through the ReLU activation function to obtain the coarse prediction sink vector.

6. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 1, characterized in that, In step S4, the Kalman filter state-space model includes a state equation and an observation equation. The state equation describes the evolution of the actual subsidence over time, and the observation equation includes a high-precision ground observation channel and a medium-precision coarse prediction observation channel.

7. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 6, characterized in that, The calculation of Kalman gain through a joint observation matrix and a piecewise filtering update strategy includes: During the period without ground-based measured data, the prior error covariance is calculated based on the state transition matrix and the posterior error covariance of the previous time step. The single-channel Kalman gain is calculated based on the prior error covariance, the coarse prediction observation matrix, and the coarse prediction observation noise covariance. The prior state estimate and the coarse prediction sinking vector are fused based on the single-channel Kalman gain to obtain the updated posterior state estimate and posterior error covariance for the single channel. During the time period in which ground-measured data exists, a joint observation matrix containing the ground observation matrix and the coarse prediction observation matrix is ​​constructed. A diagonal block joint observation noise covariance matrix containing the ground observation noise covariance and the coarse prediction observation noise covariance is also constructed. The dual-channel Kalman gain is calculated based on the prior error covariance, the joint observation matrix, and the joint observation noise covariance matrix. The prior state estimate, the ground-measured subsidence value, and the coarse prediction subsidence vector are fused based on the dual-channel Kalman gain to obtain the updated dual-channel posterior state estimate and posterior error covariance.

8. The method for predicting coal mine surrounding rock deformation based on LSTM according to claim 7, characterized in that, The process of obtaining the optimal state estimate and calculating the deviation between the optimal state estimate and the coarse predicted subsidence vector as a calibration correction includes: The posterior state estimate updated by the single channel and the posterior state estimate updated by the dual channel are used as the optimal state estimate; The difference between the optimal state estimate and the coarse predicted subsidence vector at the corresponding time is calculated to obtain the calibration correction amount; Record the calibration correction values ​​of each measuring point within a continuously set time window, construct a spatiotemporal error transfer matrix, perform principal component analysis on the spatiotemporal error transfer matrix, extract the principal components whose cumulative variance contribution rate reaches a set threshold, and obtain the spatial correlation pattern of error transfer.

Images