Surface collapse prediction method and surface collapse prediction system

Abstract

The invention discloses a surface collapse prediction method and a surface collapse prediction system, and the method comprises the steps: carrying out the synchronous monitoring of a monitoring region through employing an unmanned plane loaded with a camera and a GNSS static reference station, and obtaining the three-dimensional data of the surface of the monitoring region; establishing an oblique photography model according to the earth surface three-dimensional data; based on an oblique photography model, a prediction model is established in combination with a probability integral method and a recurrent neural network, and the prediction model at least comprises an acquisition layer, a maximum settlement height output layer and a collapse area output layer; at least predicting the maximum subsidence height and subsidence area of the surface subsidence of the monitoring area through a prediction model; therefore, the problem of manual field exploration can be solved, the problem that the image control points are difficult to arrange due to collapse is also solved, and the cost is reduced; meanwhile, high-precision prediction of the maximum ground surface settlement height and the settlement collapse area in a natural caving method mining mode is achieved through the prediction model.

Description

Technical Field

[0001] The present invention relates to the technical field of mine geological monitoring, and in particular to a surface subsidence prediction method and a surface subsidence prediction system. Background Art

[0002] Natural caving has become a common mining method when mining low-grade thick and large deposits, but this mining method can cause surface collapse, induce geological disasters, and seriously affect workers' production activities. Therefore, it is very important to monitor the mine surface in time and space continuously and predict the surface collapse.

[0003] Surface collapse prediction usually starts with establishing long-period surface time series data, among which the maximum settlement height and collapse area are important indicators. The maximum settlement height reflects the maximum vertical change during the collapse of the ore body, and the collapse area represents the settlement degree of the surface at different stages. At present, the mainstream monitoring methods include INSAR, LIDAR, GNSS and drone aerial survey methods, which can be used to extract surface settlement indicator information. However, during the monitoring process, the instability of geological structures and the hysteresis of surface settlement make it impossible to enter the collapse area, resulting in a random data observation cycle, which in turn affects the accuracy of the data. Especially when using drone aerial survey methods, a large number of ground image control points need to be deployed, and the deployment range needs to evenly cover the monitoring area, and the points are clear. However, when the surface of the mining area collapses, due to the rapid collapse rate, the deployed image control points cannot be reused. Some image control points may be damaged due to the expansion of the collapse area, which makes it difficult to continue observation. In the subsequent prediction of surface collapse, not only the prediction accuracy is affected, but also the image control points need to be re-deployed, which increases the cost. In addition, these methods also have errors caused by considering complex environments (seasons and lighting, etc.), which leads to insufficient data matching accuracy.

[0004] At the same time, there are three main methods for predicting surface subsidence: mechanical analysis, influence function method, and statistical method based on time series data. The mechanical analysis method obtains a series of complex physical parameters such as lithological characteristics and stress conditions, performs numerical simulation prediction, and then explains the mechanism of surface subsidence, but there is a certain gap between the predicted value obtained by this method and the actual value. The influence function method uses mathematical probability to explain real engineering geological events, such as probability integral method, Knothe function, etc. However, this method depends on the geological structure of the study area and is not easy to promote. The statistical method based on time series data collects a large amount of data for statistical analysis, such as random forest, support vector machine, convolutional neural network, long short-term memory network, etc. However, this method does not consider the physical mechanism of surface subsidence and lacks persuasiveness. At the same time, with the development of data analysis methods, data diversification leads to data heterogeneity, which is difficult to handle in time series prediction. Summary of the invention

[0005] The main purpose of the invention is to provide a high-precision, low-cost ground subsidence prediction method and ground subsidence prediction system.

[0006] To achieve the above object, the present invention proposes a method for predicting ground subsidence, which comprises the following steps:

[0007] Step S1, using a drone with a camera and a GNSS static reference station to synchronously monitor the monitoring area to obtain three-dimensional surface data of the monitoring area;

[0008] Step S2, establishing an oblique photography model according to the three-dimensional surface data;

[0009] Step S4: Based on the oblique photography model, a prediction model is established in combination with a probability integral method and a recurrent neural network, wherein the prediction model at least includes a collection layer, a maximum settlement height output layer, and a collapse area output layer;

[0010] Step S5: using the prediction model, at least predict the maximum settlement height and the collapse area of ​​the surface collapse in the monitoring area.

[0011] Furthermore, the step S1 specifically includes:

[0012] Step S11, using a drone with a camera to perform aerial survey over the monitoring area to obtain drone POS data, and installing the GNSS static reference station within a certain range outside the monitoring area to synchronously monitor the monitoring area to obtain ground point coordinate data;

[0013] Step S12: Based on the ground point coordinate data, the UAV POS data is processed using the PPK post-difference method to obtain the three-dimensional surface data.

[0014] Furthermore, the data obtained by the UAV and the GNSS static reference station performing synchronous monitoring of the monitoring area include at least monitoring data of 2 epochs or monitoring data of 4 satellites.

[0015] Furthermore, the step S2 specifically includes:

[0016] Step S121, calculating and obtaining the integer ambiguity;

[0017] Step S122, performing differential linear combination of the drone POS data and the ground point coordinate data to obtain virtual carrier phase observation data;

[0018] Step S123: Calculate the three-dimensional surface data according to the ground point coordinate data and the virtual carrier phase observation data.

[0019] Furthermore, after step S2 and before step S4, the following steps are further included:

[0020] Step S3: Evaluate the accuracy of the oblique photography model.

[0021] Furthermore, the step S3 specifically includes:

[0022] Step S31, selecting n checkpoints in the oblique photography model, and obtaining monitoring data of each of the checkpoints;

[0023] Step S32: Obtain the error value of the oblique photography model according to an error calculation formula, wherein the error calculation formula is: In the formula, is the measurement value of the i-th checkpoint obtained based on the oblique photography model or the drone aerial survey image, D i is the field measurement value of the i-th checkpoint;

[0024] Step S33: If the error value is within a preset range, the oblique photography model is valid; if the error value is outside the preset range, the oblique photography model is invalid, and the oblique photography model is re-established.

[0025] Furthermore, the step S4 specifically includes:

[0026] Step S41, constructing the acquisition layer based on the oblique photography model;

[0027] Step S42: Based on the acquisition layer, the maximum settlement height output layer is constructed in combination with the graph model and the RNN network;

[0028] Step S43: Based on the acquisition layer and the maximum settlement height output layer, the collapsed area output layer is constructed in combination with the SE attention mechanism and the LSTM network.

[0029] Furthermore, the step S5 specifically includes:

[0030] Step S51, obtaining the release height data of each ore pillar in the monitoring area from the oblique photography model through the acquisition layer of the prediction model;

[0031] Step S52, inputting the release height data into the maximum settlement height output layer of the prediction model for fitting, obtaining the collapse height data of each checkpoint, and obtaining the maximum settlement height data through the RNN network processing;

[0032] Step S53: input the release height data and the collapse height data into the collapse area output layer of the prediction model to calculate the collapse area of ​​the monitoring area.

[0033] The present invention also provides a surface collapse prediction system, which is applicable to the surface collapse prediction method. The surface collapse prediction system comprises:

[0034] Drones are used to carry cameras to monitor the monitoring area;

[0035] a GNSS static reference station, used to synchronously monitor the monitoring area with the UAV; and

[0036] The processing center is electrically connected to the UAV and the GNSS static reference station, and includes an oblique photography model construction module and a prediction model construction module. The oblique photography model construction module is used to establish an oblique photography model according to the three-dimensional surface data obtained by synchronous monitoring of the UAV and the GNSS static reference station. The prediction model construction module is used to establish a prediction model based on the oblique photography model in combination with a probability integral method and a recurrent neural network, and based on the prediction model, at least predict the maximum settlement height and collapse area of ​​the surface collapse in the monitoring area.

[0037] The technical solution of the present invention provides a surface collapse prediction method and a surface collapse prediction system. The surface collapse prediction method adopts a UAV remote low-altitude remote sensing method and uses a PPK post-difference method to collect monitoring area data, which can not only solve the problem of manual field exploration, but also solve the problem of difficulty in arranging image control points due to collapse, thereby reducing costs. At the same time, a probability integral method is used to fuse a recurrent neural network to construct a prediction model, and the maximum surface settlement height and settlement collapse area are predicted under the natural caving mining method. The influence of special covering layer composition on surface settlement characteristics is reduced, a multivariate RNN network is improved, the influence of internal and external factors on surface collapse is balanced, and the prediction accuracy is improved. The prediction accuracy of the maximum surface settlement height can reach 91.42%, and the prediction accuracy of the settlement collapse area can reach 87.52%. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying creative work.

[0039] Figure 1 A flow chart of an embodiment of a method for predicting ground subsidence provided by the present invention;

[0040] Figure 2 for Figure 1 Flow chart of step S1 in FIG.

[0041] Figure 3 for Figure 2 Flow chart of step S12;

[0042] Figure 4 for Figure 1 Flow chart of step S3;

[0043] Figure 5 for Figure 1 Flow chart of step S4;

[0044] Figure 6 for Figure 1 Flow chart of step S5;

[0045] Figure 7 A three-cycle oblique photography model constructed in a specific embodiment of the present invention;

[0046] Figure 8 A schematic diagram of 12 checkpoints selected in a specific embodiment of the present invention;

[0047] Fig. 9 for Figure 8 Each inspection brings you a schematic diagram of the error;

[0048] Fig.10 It is a release column model of four time periods obtained from the oblique photography model in a specific embodiment of the present invention;

[0049] Fig.11 A subsidence cloud map of a monitoring area in a time period predicted by a probability integral method in a specific embodiment of the present invention;

[0050] Fig.12 A subsidence cloud map of the monitoring area in another time period predicted by the probability integral method in a specific embodiment of the present invention;

[0051] Fig.13 A subsidence cloud map of the monitoring area in another time period predicted by the probability integral method in a specific embodiment of the present invention;

[0052] Fig.14 It is a time series diagram obtained by comparing and analyzing the predicted settlement value of the checkpoint, the PIM simulated settlement value and the actual value in a specific embodiment of the present invention;

[0053] Fig.15 A comparative analysis diagram of the predicted values ​​of the maximum periodic settlement height and the periodic collapse area and the PIM predicted values ​​in a specific embodiment of the present invention;

[0054] Fig.16 A data graph of the maximum settlement height of each cycle within a certain time period predicted by a prediction model in a specific embodiment of the present invention;

[0055] Fig.17A data graph of the collapse area of ​​each cycle predicted by a prediction model in a specific embodiment of the present invention within a certain period of time;

[0056] Fig.18 A schematic diagram of a process for predicting ground subsidence in a specific embodiment of the present invention;

[0057] Fig.19 A schematic structural diagram of an embodiment of a surface subsidence prediction system provided by the present invention.

[0058] Description of Figure Numbers:

[0059] Label name Label name 100 Surface collapse prediction system 3 Processing Center 1 Drones 31 Oblique photography model building module 2 GNSS static reference station 32 Predictive model building blocks

[0060] The realization of the purpose, functional features and advantages of the present invention will be further explained in conjunction with embodiments and with reference to the accompanying drawings. DETAILED DESCRIPTION

[0061] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work are within the scope of protection of the present invention.

[0062] It should be noted that if the embodiments of the present invention involve directional indications (such as up, down, left, right, front, back, etc.), the directional indications are only used to explain the relative position relationship, movement status, etc. between the components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication will also change accordingly.

[0063] In addition, if there are descriptions involving "first", "second", etc. in the embodiments of the present invention, the descriptions of "first", "second", etc. are only used for descriptive purposes and cannot be understood as indicating or implying their relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In addition, the meaning of "and / or" appearing in the full text includes three parallel schemes. Taking "A and / or B" as an example, it includes scheme A, or scheme B, or a scheme that satisfies both A and B. In addition, the technical solutions between the various embodiments can be combined with each other, but it must be based on the ability of ordinary technicians in the field to implement. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0064] Natural caving has become a common mining method when mining low-grade thick and large deposits, but this mining method can cause surface collapse, induce geological disasters, and seriously affect workers' production activities. Therefore, it is very important to monitor the mine surface in time and space continuously and predict the surface collapse.

[0065] Surface collapse prediction usually starts with establishing long-period surface time series data, among which the maximum settlement height and collapse area are important indicators. The maximum settlement height reflects the maximum vertical change during the collapse of the ore body, and the collapse area represents the settlement degree of the surface at different stages. At present, the mainstream monitoring methods include INSAR, LIDAR, GNSS and drone aerial survey methods, which can be used to extract surface settlement indicator information. However, during the monitoring process, the instability of geological structures and the hysteresis of surface settlement make it impossible to enter the collapse area, resulting in a random data observation cycle, which in turn affects the accuracy of the data. Especially when using drone aerial survey methods, a large number of ground image control points need to be deployed, and the deployment range needs to evenly cover the monitoring area, and the points are clear. However, when the surface of the mining area collapses, due to the rapid collapse rate, the deployed image control points cannot be reused. Some image control points may be damaged due to the expansion of the collapse area, which makes it difficult to continue observation. In the subsequent prediction of surface collapse, not only the prediction accuracy is affected, but also the image control points need to be re-deployed, which increases the cost. In addition, these methods also have errors caused by considering complex environments (seasons and lighting, etc.), which leads to insufficient data matching accuracy.

[0066] At the same time, there are three main methods for predicting surface subsidence: mechanical analysis, influence function method, and statistical method based on time series data. The mechanical analysis method obtains a series of complex physical parameters such as lithological characteristics and stress conditions, performs numerical simulation prediction, and then explains the mechanism of surface subsidence, but there is a certain gap between the predicted value obtained by this method and the actual value. The influence function method uses mathematical probability to explain real engineering geological events, such as probability integral method, Knothe function, etc. However, this method depends on the geological structure of the study area and is not easy to promote. The statistical method based on time series data collects a large amount of data for statistical analysis, such as random forest, support vector machine, convolutional neural network, long short-term memory network, etc. However, this method does not consider the physical mechanism of surface subsidence and lacks persuasiveness. At the same time, with the development of data analysis methods, data diversification leads to data heterogeneity, which is difficult to handle in time series prediction.

[0067] In view of this, the present invention provides a method for predicting surface collapse. Figures 1 to 6 A flow chart of the surface subsidence prediction method provided by the present invention.

[0068] See also Figure 1 and Fig.18, the surface collapse prediction method comprises the following steps:

[0069] Step S1: Use a drone loaded with a camera and a GNSS static reference station to synchronously monitor the monitoring area to obtain three-dimensional surface data of the monitoring area.

[0070] In this step, the drone is equipped with a camera and flies over the monitoring area according to a predetermined route, so as to obtain the original image files, camera files and POS files of the surface in the monitoring area. In addition, before conducting the aerial survey, an open area is first found within a certain range outside the monitoring area to build the GNSS static reference station. In this way, when the drone conducts aerial survey, at least one ground reference receiver in the GNSS static reference station is used to perform synchronous observation, receive real-time satellite signals, and obtain ground point coordinate data. In other words, the drone combined with the GNSS static reference station is equivalent to a mobile receiver performing real-time carrier phase observation of the satellite. In this way, the measurement error can be controlled within -6.25cm to 5.93cm, which fully meets the engineering specifications, gets rid of the constraints of image control points, overcomes the problem of inaccessibility of the monitoring area, and can replace manual measurement methods.

[0071] It should be noted that in a specific embodiment of the present invention, the monitoring area is the first mining area of ​​the Pulang Copper Mine in Yunnan. The Pulang Copper Mine in Yunnan is currently the largest mine in China that uses natural caving. It is located in the northeast of Shangri-La County, Diqing Tibetan Autonomous Prefecture in the northwest of Yunnan Province. The mining area is 72km away from the urban area of ​​Shangri-La and is under the jurisdiction of Geza Township, Shangri-La County. The mining area can be divided into a mining field, a ore dressing field and a tailings pond. The mining area is 3520-4569.6m above sea level. The first mining area of ​​the Pulang Copper Mine in Yunnan has a steep terrain, and the collapse pits formed by mining are prone to disasters such as mudslides, and manual monitoring has great safety hazards.

[0072] For further information, see Figure 2 , the step S1 specifically includes:

[0073] Step S11: Use a drone with a camera to conduct aerial survey over the monitoring area to obtain drone POS data, and install the GNSS static reference station within a certain range outside the monitoring area to synchronously monitor the monitoring area to obtain ground point coordinate data.

[0074] It should be noted that the observation equation of the GNSS static reference station for a satellite i at a certain moment is:

[0075]

[0076] In formula (1), is the carrier phase observation value of satellite i by the GNSS static reference station;

[0077] f is the carrier frequency;

[0078] c is the speed of light;

[0079] (X i , Y i , Z i ) is the spatial position of satellite i at that moment;

[0080] (X a , Y a , Z a ) is the spatial position of the GNSS static reference station at that moment;

[0081] is the ionospheric delay of the path from the GNSS static reference station to satellite i;

[0082] is the tropospheric delay of the GNSS static reference station for the satellite i path;

[0083] V ta is the clock error of the GNSS static reference station;

[0084] V ti is the clock error of satellite i;

[0085] is the integer ambiguity of the carrier phase observation of satellite i by the GNSS static reference station.

[0086] At the same time, the observation equation of the UAV to satellite i is:

[0087]

[0088] In formula (1), is the carrier phase observation value of the UAV for satellite i;

[0089] (X b , Y b , Z b ) is the spatial position of the UAV at that moment;

[0090] is the ionospheric delay of the path from UAV to satellite i;

[0091] is the tropospheric delay of the UAV to the satellite i path;

[0092] V tb is the drone’s clock error;

[0093] is the integer ambiguity of the UAV's carrier phase observation of satellite i.

[0094] Furthermore, the single difference observation equation between the GNSS static reference station and the UAV for satellite i is:

[0095]

[0096] Similarly, at the same time, the single difference observation equation between the GNSS static reference station and the UAV to another satellite j is:

[0097]

[0098] Furthermore, by subtracting equation (3) from equation (4) and performing the difference between satellite i and satellite j, the double difference observation equation of the GNSS static base station and the UAV for satellite i and satellite j can be obtained as follows:

[0099]

[0100] Based on equation (5), the GNSS static base station and the UAV simultaneously observe 5 satellites in 2 epochs, and obtain 10 double-difference observation equations, which is greater than the number of unknowns in the equation. Therefore, in order to perform PPK post-check decomposition on the UAV POS data, synchronous observation data of at least 2 epochs or more than 4 satellites is required.

[0101] That is to say, in this step, the data obtained by the UAV and the GNSS static reference station performing synchronous monitoring of the monitoring area include monitoring data of at least 2 epochs or monitoring data of at least 4 satellites.

[0102] More specifically, in a specific embodiment of the present invention, a Pegasus multi-rotor E2000S drone is used for monitoring, with an average flight speed of 14m / s, a fixed-focus single lens of 25mm, and no image control points. The flight path extends a certain distance from the ranging line, about 50 meters, to ensure that there is no leakage flight phenomenon. According to the measurement area in the north-south direction, the flight path is parallel to the map contour lines and laid in a straight line along the east-west direction. The take-off point is set near the office area at an altitude of 3720m. A reference station is set next to the take-off point for synchronous GNSS static observation during aerial measurement. According to the highest altitude of the surface subsidence area of ​​4100 meters, the take-off point is 1.5 kilometers away from the farthest point of the survey area, the flight altitude is set to 500 meters, the total flight time is 25 minutes, and there are a total of two sorties.

[0103] In addition, the GNSS static reference station is built within a range of 2 kilometers outside the monitoring area, and the UAV and the GNSS static reference station perform synchronous monitoring of the monitoring area.

[0104] Step S12: Based on the ground point coordinate data, the UAV POS data is processed using the PPK post-difference method to obtain the three-dimensional surface data.

[0105] For further information, see Figure 3 , the step S12 specifically includes:

[0106] Step S121, calculate and obtain the integer ambiguity.

[0107] In this step, the integer ambiguity is calculated by initializing the OTF (Optical Transfer Function).

[0108] Step S122: performing a differential linear combination of the drone POS data and the ground point coordinate data to obtain virtual carrier phase observation data.

[0109] In this step, the relative position relationship between the GNSS static reference station and the UAV is obtained through the virtual carrier phase observation data.

[0110] Step S123: Calculate the three-dimensional surface data according to the ground point coordinate data and the virtual carrier phase observation data.

[0111] Step S2: establishing an oblique photography model based on the three-dimensional surface data.

[0112] In this step, a model is built based on the three-dimensional surface data through modeling software to obtain an oblique photography model (real-scene three-dimensional model) for restoring the surface appearance.

[0113] Specifically, the modeling software includes Street Factory from France's Air Bus, PhotoScan from Russia's Agisoft, Pix4D mapper from Switzerland's Pix4D, Context Capture from Bentley, and domestic DJI Intelligent Mapping.

[0114] More specifically, in the specific embodiments of the present invention, please refer to Figure 7 , the complete surface subsidence information and related mining data for 35 months from February 2018 to January 2021 were collected, and the oblique photography model was established with a monthly cycle, where a, b, and c represent the surface information on January 17, 2019, January 17, 2020, and January 17, 2021, respectively.

[0115] Step S4: Based on the oblique photography model, a prediction model is established in combination with a probability integral method and a recurrent neural network, wherein the neural network of the prediction model includes at least a collection layer, a maximum settlement height output layer, and a collapse area output layer.

[0116] Step S5: using the prediction model, at least predict the maximum settlement height and the collapse area of ​​the surface collapse in the monitoring area.

[0117] For details, please refer to Figure 1 , after step S2 and before step S4, further comprising:

[0118] Step S3: Evaluate the accuracy of the oblique photography model.

[0119] In this step, the prediction accuracy of the subsequent prediction model is improved by determining the accuracy of the oblique photography model.

[0120] For further information, see Figure 4 , the step S3 specifically includes:

[0121] Step S31, selecting n checkpoints in the oblique photography model, and obtaining monitoring data of each of the checkpoints.

[0122] It should be noted that the checkpoints may also be selected through drone aerial survey images.

[0123] Step S32: Obtain the error value of the oblique photography model according to an error calculation formula, wherein the error calculation formula is:

[0124]

[0125] In formula (6), is the measurement value of the i-th checkpoint obtained based on the oblique photography model or the drone aerial survey image;

[0126] D i is the field measurement value of the i-th checkpoint.

[0127] It should be noted that the measurement value includes the elevation value and the horizontal value of each checkpoint.

[0128] Step S33: If the error value is within a preset range, the oblique photography model is valid; if the error value is outside the preset range, the oblique photography model is invalid, and the oblique photography model is re-established.

[0129] In this step, in one embodiment of the present invention, when evaluating the accuracy of the elevation value of each checkpoint, the preset range is ±0.3m, and when evaluating the accuracy of the horizontal value of each checkpoint, the preset range is ±0.5m.

[0130] like Figure 8 As shown in the specific embodiment, 12 relatively fixed and non-collapsed checkpoints (K1-K12) are selected in the subsidence area of ​​Pulang Copper Mine in Yunnan Province, and the checkpoint error is Fig. 9 It can be seen that the errors of DX, DY, and DZ are controlled within -6.25cm-5.93cm, R x The value is 2.493cm, R yThe value is 2.732cm, R z The value is 3.448 cm, which meets the specifications of geographic information engineering (Bureau 2012).

[0131] For further information, see Figure 5 , the step S4 specifically includes:

[0132] Step S41: construct the acquisition layer based on the oblique photography model.

[0133] In this step, in the acquisition layer, the monitoring area is extracted through the oblique photography model to obtain the monitoring data of the checkpoint, wherein a part of the monitoring data forms a time series data set, including the time series data of all ore-releasing points composed of monitoring data such as the ore-releasing amount, sinking coefficient, pillar size, maximum impact angle and release height of each ore-releasing point in the monitoring time period, which is used for the subsequent maximum settlement height prediction, and also includes the number of ore-releasing points, maximum settlement height, bottom pulling area and backfilling amount of the collapse surface in the monitoring area each month, which is used for the subsequent collapse area prediction. The other part of the monitoring data is used as a verification data set to evaluate the accuracy of the subsequent prediction results.

[0134] It should be noted that the acquisition layer can also obtain monitoring data corresponding to multiple checkpoints based on the images obtained by the drone aerial survey.

[0135] Step S42: Based on the acquisition layer, the maximum settlement height output layer is constructed in combination with the graph model and the RNN network.

[0136] It should be noted that the amount of collapse in the monitoring area is greatly affected by the collapse and movement of underground ore and rock, which is manifested in that the amount of ore released is positively correlated with the collapse. Therefore, the maximum surface settlement height predicted for each cycle occurs above the center of the ore pillar released in that cycle. In this regard, the checkpoint also includes the ground point directly above the ore release point in the collapse surface observed in each cycle. By judging the settlement height of each checkpoint in each cycle, the maximum settlement height of the collapse surface in each cycle can be extracted.

[0137] In this step, the probability integral method is used to process the time series data set, that is, the ore placement height data is processed to obtain the ore placement height, which is input into multiple graph unit structures of the graph model as the feature vector of each checkpoint and enters the RNN network.

[0138] For details, please refer to Fig.18, the graph unit structure is divided into nodes and edges. The node features are updated through the adjacency matrix. During this period, the time information is added to the input embedding for position embedding. The temporal nature of the nodes is strengthened through position encoding. The edge updates the weight parameters through the self-attention mechanism between nodes. It is a channel for information transmission between connecting nodes and has directionality. A fully connected network can be connected, and mapping is performed through cross attention. The lower the attention, the more effectively suppressed, and the weights between all nodes are obtained.

[0139] A single graph unit structure can extract feature information of a single time point, use the RNN network structure to establish data features between moments, iterate through two layers of RNN networks, and strengthen the features between moments through a multi-head attention mechanism. This process updates the overall weights of each feature of the graph unit structure, and the resulting output vector is connected to a point-to-point feedforward network for one-step processing, and the output is once again connected to a residual network for normalization. This process is similar to the Encoder processing in the transformer, and its purpose is to express the input continuously with attention information, which helps to better learn the features. The output value is connected to the fully connected network and performs SoftMax regression calculations, so as to evaluate the accuracy of the model value, that is, the data in the verification data set, and the maximum value is selected as the output of this layer.

[0140] Step S43: Based on the acquisition layer and the maximum settlement height output layer, the collapsed area output layer is constructed in combination with the SE attention mechanism and the LSTM network.

[0141] In this step, the data output by the acquisition layer and the maximum settlement height output layer are processed by the SE (Squeeze-and-Excitation) attention mechanism, and then the processed data is imported into the LSTM network for iterative calculation to obtain the collapsed area data, and then the collapsed area data is normalized by the ResNet residual network, and finally the model value, i.e., the data in the validation data set, is evaluated for accuracy and output. In this way, the problem of long sequence prediction memory loss and multi-layer regression causing prediction distortion in the LSTM network can be solved.

[0142] It should be noted that the SE (Squeeze-and-Excitation) attention mechanism is an attention module that performs attention operations on the channel dimension. The purpose of this operation is to focus on the attention weights of inter-channel features and suppress unimportant channel features. Its implementation steps are generally divided into three steps: first, the Squeeze compression operation compresses each channel feature and keeps the number of feature channels unchanged, and globally pools each channel feature to convert it into a real value; second, Excitation activation, which maps the global features of two fully connected layer networks into a weight vector and normalizes it using the sigmoid function to obtain the attention weight; finally, the attention weight is multiplied by the original feature to obtain a new feature, which can be imported into the LSTM network for time series data feature extraction.

[0143] After passing through the LSTM network, the data is superimposed with the PIM output value through transmission, and more feature information is extracted through convolution. The PIM residual network is based on the RNN output value. The network structure performs Resnet residual network normalization calculation on the RNN output value and the superimposed output value, and its purpose is to screen the feature value. This is because there are two layers of regression operations in the RNN output value, and the error itself is large. At this time, a set of collapsed area values ​​is obtained by fitting the output value of the probability integral method. This value is superimposed with the RNN output value, and the residual network calculation is introduced to obtain more feature information, reduce the error caused by multi-layer regression, and approach the true value.

[0144] In this step, the surface collapse conditions at least include the maximum settlement height and the collapse area of ​​the surface collapse.

[0145] For further information, see Figure 6 , the step S5 specifically includes:

[0146] Step S51, obtaining the release height data of each ore pillar in the monitoring area from the oblique photography model through the acquisition layer of the prediction model.

[0147] Step S52: input the release height data into the maximum settlement height output layer of the prediction model for fitting to obtain the collapse height data of each checkpoint, and obtain the maximum settlement height data through the RNN network processing.

[0148] In this step, the accuracy of the obtained collapse height data is low, and the influence of external factors such as weather, rainfall and artificial backfilling on the change of surface collapse is not considered. Therefore, it needs to be processed by the RNN network to improve the data accuracy.

[0149] Step S53: input the release height data and the collapse height data into the collapse area output layer of the prediction model to calculate the collapse area of ​​the monitoring area.

[0150] More specifically, in a specific embodiment of the present invention, the mining level of the monitoring area is 3720 meters above sea level, which is used as the horizontal working surface. According to the actual situation, the surface collapse of the Pulang Copper Mine began in the 12th month after it was put into production, that is, the surface began to collapse on February 16, 2018. Therefore, the prediction of surface collapse caused by mining starts from February 17, 2018. At this time, there are 121 active mining points, and the number of mining points gradually increases every day. The time period for surface collapse prediction is selected as one month. According to the expected collapse period, the input mining points and mining output are the actual mining points and mining output of this month, and the calculation results of this month are used as the calculation input surface of the next month.

[0151] Furthermore, through the acquisition layer of the prediction model, the model of the pillar released in this cycle is obtained from the oblique photography model, and the release height data of the pillar can be obtained through simulation, such as Fig.10 As shown, the release column model of the four selected time periods is generated, from which the change law of the pillar can be seen.

[0152] Furthermore, the release height data is input into the maximum settlement height output layer of the prediction model for fitting, and only the probability integral method (PIM) is used to predict the settlement value of each checkpoint, that is, the collapse height data of each checkpoint. Figure 11-Figure 3 As shown, it is a collapse cloud map generated by fitting, which shows the vertex sinking amount in each cycle after simulation, and draws contour lines, where the vertex sinking amount can be regarded as the collapse height data of the checkpoint. Fig.11 The monthly collapse cloud map from February 2018 to January 2019 shows that the amount of subsidence at the top slows down from the central collapse area to the surrounding areas. Fig.11 -d and Fig.11 -e The relative amount of sinking of the vertices in the central area Fig.11 -a, Fig.11 -b and Fig.11 -c surges, Fig.11 -f The amount of subsidence at the apex of the central area decreased, and the collapse range expanded to the north. Fig.11 -h to Fig.11 -k The sinking value of the vertex in the central area increases, with a value of 4-5m, and the degree of collapse gradually increases from the central area to the surrounding areas; Fig.12 It shows the monthly subsidence cloud map from January 2019 to January 2020, where Fig.12 -b and Fig.12 -I has a lower collapse degree, and the maximum peak subsidence in other months is maintained between 3.9-6.3m. Fig.12 -f shows that starting from June 2019, the maximum peak subsidence gradually shifted from the southern part of the central area to the northern part; Fig.13 It shows the monthly subsidence cloud chart from January 2019 to January 2020. Fig.13 -a to Fig.13 -c The maximum vertex sinking amount is offset to the east, Fig.13 -d does not sink much, Fig.13 -e to Fig.13 -h The maximum peak sinking is maintained at 4.8-6.3m, Fig.13 -I to Fig.13 -LThe maximum vertex subsidence shifts to the northwest, and the value gradually decreases, but the overall collapse range increases.

[0153] from Figure 11-13 The maximum settlement height value of each cycle and the settlement location distribution of the cycle can be seen in it. However, the accuracy of the prediction results using only the probability integral method is low, and the influence of external factors such as weather, rainfall and artificial backfilling on the change of surface collapse is not considered. Therefore, the results are imported into the RNN network for processing. Specifically, the data is divided into three test data, namely February 2018 to February 2019, February 2019 to March 2020, and March 2020 to January 2021. Two ore-dropping points are selected in each set, which are all newly added ore-dropping points. The settlement prediction value of the checkpoint, the PIM simulated settlement value and the actual value are compared and analyzed, and drawn into a time series diagram (such as Fig.14 and integrate the maximum settlement height data of the period, and then use the determination coefficient (R 2 ), root mean square error (RMSE) and mean absolute percentage error (MAPE) were used for evaluation, and various indicators were compared through prediction by different models.

[0154] Furthermore, the maximum settlement height data of each cycle is used as input, and the samples are the monthly surface collapse surfaces, and the collapse area of ​​each cycle is obtained. The predictions of the maximum settlement height and collapse area of ​​each cycle based on this model are compared with the PIM prediction values, such as Fig.14 and Fig.15 As shown. Fig.16 It can be seen that the PIM prediction results are generally lower than the actual settlement values, especially during the period from July to December, when the gap between the PIM prediction method and the actual value is obvious, which is consistent with the rainy season in the region. Fig.17 It can be seen that from February 2018 to April 2019, there was a huge gap between the PIM prediction method and the actual value. During this period, artificial backfilling was carried out to slow down the speed of surface collapse, and the PIM prediction method did not take into account the impact of this factor.

[0155] In the technical solution of the present invention, a UAV remote low-altitude remote sensing method is adopted, and the PPK post-difference method is used to collect monitoring area data, which can not only solve the problem of manual field exploration, but also solve the problem of difficulty in arranging image control points due to collapse, thereby reducing costs; at the same time, a probability integral method is adopted to fuse a recurrent neural network to construct a prediction model, and the maximum surface settlement height and settlement collapse area are predicted under the natural caving mining method. It reduces the influence of special covering layer composition on surface settlement characteristics, improves the multivariate RNN network, balances the influence of internal and external factors on surface collapse, and improves the prediction accuracy. The prediction accuracy of the maximum surface settlement height can reach 91.42%, and the prediction accuracy of the settlement collapse area can reach 87.52%.

[0156] The present invention also provides a surface collapse prediction system 100, which is applicable to the surface collapse prediction method described above. Fig.18 The surface subsidence prediction system includes an unmanned aerial vehicle 1, a GNSS static reference station 2 and a processing center 3. The unmanned aerial vehicle 1 is used to carry a camera to perform aerial survey of the monitoring area; the GNSS static reference station 2 is used to synchronously monitor the monitoring area with the unmanned aerial vehicle 1; the processing center is electrically connected to the unmanned aerial vehicle 1 and the GNSS static reference station 2, and includes an oblique photography model construction module 31 and a prediction model construction module 32. The oblique photography model construction module 31 is used to establish an oblique photography model based on the three-dimensional surface data obtained by synchronously monitoring the unmanned aerial vehicle 1 and the GNSS static reference station 2. The prediction model construction module 32 is used to establish a prediction model based on the oblique photography model and in combination with the probability integral method and the recurrent neural network, and based on the prediction model, at least predict the maximum settlement height and collapse area of ​​the surface subsidence in the monitoring area.

[0157] The above description is only a preferred embodiment of the present invention, and does not limit the patent scope of the present invention. All equivalent structural changes made by using the contents of the present invention specification and drawings under the inventive concept of the present invention, or directly / indirectly applied in other related technical fields are included in the patent protection scope of the present invention.

Claims

1. A method for predicting surface subsidence, characterized in that: The surface subsidence prediction method comprises the following steps: Step S1, using a drone with a camera and a GNSS static reference station to synchronously monitor the monitoring area to obtain three-dimensional surface data of the monitoring area; Step S2, establishing an oblique photography model according to the three-dimensional surface data; Step S4: Based on the oblique photography model, a prediction model is established in combination with a probability integral method and a recurrent neural network, wherein the prediction model at least includes a collection layer, a maximum settlement height output layer, and a collapse area output layer; Step S5: using the prediction model, at least predict the maximum settlement height and the collapse area of ​​the surface collapse in the monitoring area.

2. The method for predicting ground subsidence according to claim 1, characterized in that: The step S1 specifically includes: Step S11, using a drone with a camera to perform aerial survey over the monitoring area to obtain drone POS data, and installing the GNSS static reference station within a certain range outside the monitoring area to synchronously monitor the monitoring area to obtain ground point coordinate data; Step S12: Based on the ground point coordinate data, the UAV POS data is processed using the PPK post-difference method to obtain the three-dimensional surface data.

3. The method for predicting ground subsidence according to claim 2, characterized in that: The data obtained by the UAV and the GNSS static reference station performing synchronous monitoring of the monitoring area include at least monitoring data of 2 epochs or monitoring data of 4 satellites.

4. The method for predicting ground subsidence according to claim 2, characterized in that: The step S2 specifically includes: Step S121, calculating and obtaining the integer ambiguity; Step S122, performing differential linear combination of the drone POS data and the ground point coordinate data to obtain virtual carrier phase observation data; Step S123: Calculate the three-dimensional surface data according to the ground point coordinate data and the virtual carrier phase observation data.

5. The method for predicting ground subsidence according to claim 4, characterized in that: After step S2 and before step S4, the method further includes: Step S3: Evaluate the accuracy of the oblique photography model.

6. The method for predicting ground subsidence according to claim 5, characterized in that: The step S3 specifically includes: Step S31, selecting n checkpoints in the oblique photography model, and obtaining monitoring data of each of the checkpoints; Step S32: Obtain the error value of the oblique photography model according to an error calculation formula, wherein the error calculation formula is: In the formula, is the measurement value of the i-th checkpoint obtained based on the oblique photography model or the drone aerial survey image, D i is the field measurement value of the i-th checkpoint; Step S33: If the error value is within a preset range, the oblique photography model is valid; if the error value is outside the preset range, the oblique photography model is invalid, and the oblique photography model is re-established.

7. The method for predicting ground subsidence according to claim 6, characterized in that: The step S4 specifically includes: Step S41, constructing the acquisition layer based on the oblique photography model; Step S42: Based on the acquisition layer, the maximum settlement height output layer is constructed in combination with the graph model and the RNN network; Step S43: Based on the acquisition layer and the maximum settlement height output layer, the collapsed area output layer is constructed in combination with the SE attention mechanism and the LSTM network.

8. The method for predicting ground subsidence according to claim 7, characterized in that: The step S5 specifically includes: Step S51, obtaining the release height data of each ore pillar in the monitoring area from the oblique photography model through the acquisition layer of the prediction model; Step S52, inputting the release height data into the maximum settlement height output layer of the prediction model for fitting, obtaining the collapse height data of each checkpoint, and obtaining the maximum settlement height data through the RNN network processing; Step S53: input the release height data and the collapse height data into the collapse area output layer of the prediction model to calculate the collapse area of ​​the monitoring area.

9. A surface subsidence prediction system, applicable to the surface subsidence prediction method according to any one of claims 1 to 8, characterized in that: The surface subsidence prediction system comprises: Drones are used to carry cameras to monitor the monitoring area; a GNSS static reference station, used to synchronously monitor the monitoring area with the UAV; and The processing center is electrically connected to the UAV and the GNSS static reference station, and includes an oblique photography model construction module and a prediction model construction module. The oblique photography model construction module is used to establish an oblique photography model according to the three-dimensional surface data obtained by synchronous monitoring of the UAV and the GNSS static reference station. The prediction model construction module is used to establish a prediction model based on the oblique photography model in combination with a probability integral method and a recurrent neural network, and based on the prediction model, at least predict the maximum settlement height and collapse area of ​​the surface collapse in the monitoring area.

Images