A kind of open-pit to underground mining ground subsidence monitoring method based on point cloud double-phase difference

Abstract

This invention proposes a surface subsidence monitoring method for open-pit to underground mining transition based on point cloud dual-temporal difference. Belonging to the interdisciplinary field of artificial intelligence, mining, and safety, it aims to solve the technical challenges of severe information loss in 2D methods and insufficient accuracy and reliance on extensive manual annotation in 3D methods for surface subsidence monitoring during open-pit to underground mining transitions. The invention uses an SKPConv-Net neural network as its core, directly processing raw 3D point cloud data through a dual-branch encoder and decoder architecture. The encoder employs kernel-point convolution with rigid kernels to extract multi-scale features from the dual-temporal point cloud, calculates feature differences using a nearest neighbor matching strategy, and inputs these features into the decoder, ultimately outputting the point-scale subsidence change results. By pre-training the network using simulated datasets, the need for annotation of real mine data is reduced. High-precision monitoring is achieved in scenarios with multiple sensors, high noise, and point clouds of varying densities, providing efficient technical support for dynamic monitoring and safety management of surface subsidence in complex open-pit to underground mining scenarios.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of artificial intelligence, mining, and safety, and specifically relates to a method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference. Background Technology

[0002] In open-pit to underground mining operations, goaf collapse is a typical geological hazard that accompanies resource extraction. Surface subsidence is one of the important monitoring methods for goaf collapse, and its monitoring and early warning are crucial for mine safety. Currently, surface subsidence monitoring technologies are mainly divided into methods based on 2D data and traditional methods based on 3D point clouds. However, both types of methods have limitations such as information loss and insufficient accuracy. At the same time, the development of deep learning in the field of 3D point cloud processing also faces bottlenecks in adapting to mining scenarios, as detailed below: Traditional 2D monitoring methods largely rely on optical imagery and digital surface models (DSMs). They identify subsidence areas by comparing pixel differences or elevation changes across different time phases of data. For example, DSM differential calculations are used to determine elevation changes, which are then combined with threshold segmentation to define the subsidence range. However, these methods inevitably lose crucial 3D information such as building facades and surface micro-topography during the rasterization process, making it impossible to accurately depict the spatial distribution characteristics of subsidence in complex mine terrain. Furthermore, 2D data is prone to spurious changes due to factors such as lighting conditions, weather variations, and sensor viewing angle differences. This is especially true in environments with high dust levels and unstable lighting, further reducing monitoring accuracy. While traditional methods based on 3D point clouds can avoid the information loss problem of 2D methods, after acquiring the original 3D point cloud through devices such as LiDAR, the ICP algorithm is used to register the point cloud and calculate the distance between point clouds to identify settlement points. Alternatively, settlement areas can be classified based on manually designed features such as normal vectors and curvature, combined with traditional machine learning algorithms such as random forests. However, these methods are highly dependent on manually designed features, have poor generalization ability, and are difficult to adapt to complex scenarios in mines with large differences in point cloud density and high noise levels. Furthermore, the sparsity and disorder of point clouds lead to low registration accuracy, ultimately resulting in large settlement classification errors, which cannot meet the high-precision monitoring needs of mines.

[0003] In recent years, deep learning technology has made significant progress in areas such as 3D point cloud semantic segmentation and target monitoring. However, in the transition from open-pit to underground mining, surface subsidence is affected by both open-pit and underground mining, posing multiple challenges to its monitoring. First, existing 3D point cloud deep learning networks are mostly designed for general scenarios such as urban building segmentation and autonomous driving target monitoring, without considering the special characteristics of mine point clouds. Mine point clouds have a single target type and fluctuate greatly in data quality. Directly migrating general networks to the mine scenario will result in insufficient feature extraction and a significant decrease in monitoring accuracy. Second, deep learning methods typically require large-scale labeled bi-temporal point cloud data to classify each point as having no subsidence, slight subsidence, moderate subsidence, or severe subsidence. However, labeling real mine point clouds requires professionals to combine geological exploration reports and GNSS field measurement data, which is not only time-consuming and labor-intensive, but the labeling cost for a single scenario can even reach hundreds of thousands of yuan. It is difficult to form a large-scale labeled dataset that meets the training requirements, which seriously limits the practical application of deep learning methods in mines. In addition, mine data often comes from different sensors such as satellite photogrammetry, UAV LiDAR, and ground 3D laser scanning. Existing networks mostly adopt a shared weight design, which cannot effectively adapt to the feature differences of multi-sensor data. In cross-sensor data fusion monitoring scenarios, the problem of sharp drop in accuracy is particularly prominent, which further restricts the promotion and application of the technology.

[0004] In summary, existing surface subsidence monitoring technologies for open-pit to underground mining generally suffer from problems such as insufficient information integrity, weak generalization ability, high annotation costs, and poor adaptability to multiple sensors. There is an urgent need for a high-precision monitoring method that can directly process raw 3D point clouds, adapt to complex monitoring scenarios in open-pit to underground mining, and reduce annotation dependence, so as to meet the actual needs of dynamic monitoring and safety management of surface subsidence in open-pit mining sites. Summary of the Invention

[0005] The purpose of this invention is to provide a method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference, which solves the problems of information loss in existing 2D methods, weak generalization ability of traditional 3D methods, high annotation cost of deep learning methods, and poor adaptability to multiple sensors, so as to achieve accurate, efficient and low-cost monitoring of surface subsidence in open-pit to underground mining.

[0006] The core technical solution of this invention revolves around the construction of the SKPConv-Net neural network, the design of a mining-specific point cloud dataset, and model training optimization strategies. Specifically, it includes three key steps: First, designing an SKPConv-Net neural network that integrates a dual-branch symmetric structure with a rigid kernel KPConv, constructing a dual-branch encoder and decoder architecture to achieve end-to-end settlement change monitoring of the original 3D point cloud; Second, constructing a point cloud dataset containing both simulated and real datasets for the open-pit to underground mining scenario, providing data support for model training and optimization; Third, adopting a training strategy of pre-training with simulated data plus fine-tuning with real data to reduce the dependence on real data annotation, while adapting to multi-sensor data scenarios through two network operation modes, ultimately achieving the classification of settlement changes in the dual-temporal point cloud to be monitored.

[0007] To achieve the above objectives, the technical solution adopted in this invention is: a method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference, the steps of which are as follows: Step 1: Preprocessing of dual-temporal point cloud data for outdoor to underground scenes: Data preprocessing is the foundation for ensuring effective feature extraction by SKPConv-Net. Addressing issues such as dust interference in outdoor to underground point cloud data, inconsistencies in coordinates due to sensor differences, and uneven point cloud density, the following preprocessing steps are designed: 1.1 Noise Removal: Input dual-temporal original 3D point cloud (Time t1, historical data) and (At time t2, current data) A statistical filtering algorithm is used to remove outlier noise points. Specifically, for each point, the average distance of its K nearest neighbors (K=20) is calculated. If the deviation of the point from the average distance exceeds a preset standard deviation threshold (threshold=1.0), it is identified as a noise point and removed. The formula is as follows:

[0008] in, For point To its The average distance between nearest neighbors, The mean of the average distances between all points. The standard deviation of the mean distance This is the standard deviation threshold.

[0009] 1.2 Coordinate System 1: Since the dual-temporal point clouds may originate from different sensors or different acquisition periods, resulting in coordinate reference differences, the Iterative Closest Point (ICP) algorithm is used to register the dual-temporal point clouds to the same coordinate system. The origin is the mine's pre-set benchmark control point, with the X-axis pointing in the main mining direction and the Z-axis representing the elevation direction. The ICP algorithm solves for the optimal rigid transformation matrix by minimizing the Euclidean distance error between corresponding points in the dual-temporal point clouds. ( For rotation matrix, (where the vector is a translation vector), the formula is as follows:

[0010] in, , Corresponding point pairs in the two-phase point cloud are determined through nearest neighbor search. After registration, Convert to Ensure with The spatial locations are consistent.

[0011] 1.3 Input Feature Construction: If the original point cloud contains supplementary information such as LiDAR intensity and RGB, it is fused with 3D coordinates (X,Y,Z) as the input feature of this network to enhance feature discriminativeness; if it only contains coordinate information, the 3D coordinates are directly used as the input feature, and the input feature dimension is reduced. .

[0012] 1.4 Sub-point cloud segmentation: Considering the large scale of point clouds in open-pit mines, directly inputting them into the network would lead to insufficient computing resources. Therefore, the registered dual-temporal point clouds are segmented... and The system divides the data into vertical cylindrical sub-regions, with the cylinder radius set according to the density of the mine's point cloud: 30m for sparse data and 25m for dense data. This ensures that each sub-point cloud contains sufficient surface detail and potential subsidence areas, while avoiding computational inefficiencies caused by excessively large sub-point clouds. Adjacent sub-point clouds maintain a 30% overlap during segmentation to prevent subsidence areas from being truncated and to ensure the continuity of subsequent monitoring results.

[0013] Step 2: SKPConv-Net Neural Network Design: A dual-branch encoder and decoder architecture is adopted. Point cloud features are extracted based on kernel-level convolution (KPConv), and feature difference calculation and skip connections are combined to achieve end-to-end settlement change monitoring. The overall network structure includes an input layer, a dual-branch encoder, a feature difference calculation layer, a decoder, and an output layer. The design of each module is as follows: 2.1 Input Layer Module: Receives the preprocessed dual-temporal sub-point cloud features and coordinates from 1); for each segmented sub-point cloud pair... , (k is the sub-point cloud index), the input feature is and ( , (These are the number of points in the two sub-point clouds, respectively), and their corresponding coordinates are... and The data can be directly input into the dual-branch encoder without the need for rasterization or voxelization preprocessing, thus avoiding information loss.

[0014] 2.2 Dual-branch encoder module: The encoder adopts a dual-branch structure, with the two branches having completely identical architectures, each processing... and Each branch contains 5 layers, each consisting of 2 convolutional blocks. Except for the first layer, the first convolutional block of the remaining layers employs a kernel-based convolutional downsampling design with stride. Multi-scale feature extraction is achieved through recursive downsampling, while maintaining robustness to point clouds of varying densities in the mine. KPConv captures local geometric features by defining kernel distribution and correlation functions within the point cloud neighborhood. Its core formula includes two parts: convolution output calculation and kernel function definition. For point clouds... For any point x in the matrix, the convolution result is the weighted sum of the features of all points in the neighborhood of that point, as shown in the formula:

[0015] in It is a three-dimensional neighborhood centered at x with radius R. Neighboring points Input features, For kernel functions; kernel function g is constructed based on preset kernel points and related functions, and it is applied to neighboring points. offset vector relative to center x Its form is (K is the number of core points, To preset the spatial location of the core point, (for learnable weight matrix), correlation function The distance affected by the core point is positively correlated with R, thus ensuring stable feature extraction in sparse point cloud scenarios.

[0016] To achieve multi-scale feature extraction, the encoder performs recursive downsampling through kernel convolution with stride, with each layer's cell size... With convolution radius according to Recursive definition ( (The initial cell size is set based on the point cloud density). Through this design, the encoder progressively reduces the number of points in the point cloud (downsampling rate of 2) and increases the feature dimension (from...) from layer 1 to layer 5. Up to 512 dimensions), the final output is a multi-scale feature matrix of the dual-temporal sub-point cloud: Branch 1 processing Output and Branch 2 processing Output and Meanwhile, the intermediate features of layers 1-4 of the encoder are passed to the decoder through skip connections to supplement the details lost during downsampling.

[0017] 2.3 Feature Difference Calculation Layer Module: Even in areas without settlement, the number of points in a dual-temporal point cloud can mismatch with their spatial location due to differences in sensor scanning angle and point cloud density, making direct calculation of feature differences prone to errors. Therefore, this layer adopts a nearest neighbor matching-based strategy to calculate the differences in the dual-temporal feature matrices. The specific steps are as follows: right Each point in ,exist Find the point with the closest spatial distance. ,Right now ;extract corresponding feature vector and corresponding feature vector The difference between the two is calculated as a characteristic of change, using the following formula:

[0018] in, The feature difference matrix centrally reflects the settlement-related changes in the two-phase point cloud and serves as the core basis for subsequent settlement classification.

[0019] 2.4 Rotation or Deformation Analysis Module

[0020] To accurately quantify the rotation and non-rigid deformation characteristics of dual-temporal point clouds and provide supplementary geometric basis for settlement classification, a rotation or deformation analysis module has been added after feature difference calculation. This module, based on multi-scale features extracted by KPConv and registered coordinate information, completes the quantitative calculation of rotation angle and deformation intensity in two steps.

[0021] First, the rotation parameters are solved, specifically for the feature difference matrix. The selected candidate regions of change are analyzed based on the intensity of feature differences. Sure, Let be the 90th quantile of the feature differences in the dataset, assuming that the displacement of points within the region conforms to the law of rigid body rotation. For each point within the candidate region... Reuse the kNN matching results from feature difference calculation to obtain its... nearest neighbor in Establish the coordinate correspondence. Solve the rotation matrix R and translation vector t using least squares optimization, with the optimization objective being:

[0022] Where M is the number of points in the candidate region, It is an orthogonal rotation matrix. This is a translation vector. To solve, first calculate the coordinates of the center point of the cloud at the two time points. , After eliminating the translation effect, construct the covariance matrix. ,right Perform SVD decomposition Then the rotation matrix is Rotation matrix Convert to Euler angles ( These correspond to rotation angles around the X, Y, and Z axes, respectively, to quantify the overall rotation trend of the candidate region. The formula is as follows:

[0023] Next, deformation strength calculation is performed. For non-rigid deformation regions, a comprehensive evaluation is conducted using normal vector offset and local strain. The normal vectors at corresponding points at times t1 and t2 are calculated. , Define the normal vector offset Simultaneously, a local neighborhood strain tensor is constructed for each point. Given the set of k nearest neighbors, calculate the covariance matrix of the displacement vector. The principal strain is obtained through eigenvalue decomposition. Define the deformation intensity Ultimately, the rotation angle ( ) and deformation intensity Concatenate into geometric feature vectors , and the feature difference matrix The fusion yields an enhanced feature matrix. This provides more comprehensive feature support for subsequent decoder classification.

[0024] 2.5 Decoder Module: The decoder adopts a 5-layer network structure symmetrical to the encoder. Through a combination of nearest neighbor upsampling, skip connection fusion, and unary convolution, it gradually restores the point cloud resolution while fusing intermediate features from the encoder, ensuring that the output point-scale sinking classification results retain sufficient detail. The specific design is as follows: First, nearest neighbor upsampling is performed on the feature difference matrix. Perform upsampling, setting the sampling rate to 2 for each layer, and use nearest neighbor interpolation to adjust the feature dimensions from the downsampled values. Restore to the original sub-point cloud point count First, it ensures that the final classification result corresponds one-to-one with the original point cloud; second, it performs skip connection fusion. After each layer completes upsampling, the intermediate features of the corresponding layer of the encoder are concatenated with the upsampled features to supplement the local details lost during downsampling. The formula for this process is as follows:

[0025] in Represents the features of the j-th layer of the decoder. Represents the encoder's first Intermediate features of the layer For nearest neighbor upsampling operation, The first step involves feature concatenation; the final step is univariate convolution and nonlinear transformation. For each layer of fused features, univariate convolution is used for dimensionality compression and feature transformation, gradually mapping the feature dimension from 512 dimensions to the settlement category dimension (no settlement, slight settlement, moderate settlement, severe settlement). Simultaneously, the ReLU activation function is introduced to enhance the nonlinear expressive power of the features. The formula is as follows:

[0026] in , These are the weight matrix and bias vector of the j-th unary convolution, respectively.

[0027] 2.6 Output Layer Module: The output layer performs softmax normalization on the transform features of the last layer of the decoder to obtain the sedimentation class probability for each point, as shown in the following formula:

[0028] in, The i-th point belongs to the c-th type of settlement ( No settlement. Slight settlement Moderate settlement. The probability of severe subsidence. Let be the feature value of the i-th point and c-th dimension in the last layer of the decoder.

[0029] To address the problem of class imbalance in mine settlement data, a point-level negative log-likelihood loss function (NLL) is used to optimize the network parameters, as shown in the following formula:

[0030] in, Let i be the actual settlement category label for the i-th point. Here, N represents the predicted probability for the corresponding category, and N is the number of points in the sub-point cloud.

[0031] 2.7 Network Operation Modes: To adapt to multi-sensor data from open-pit to underground mining, SKPConv-Net is designed with two operation modes: Shared weight mode: The dual-branch encoder shares the same network weights and hyperparameters, which is suitable for dual-temporal point clouds with the same sensor and density. This mode has fewer parameters and only needs to train one set of encoder weights, resulting in high training efficiency. It can make full use of the feature consistency of data from the same sensor to improve monitoring accuracy. Non-shared weight mode: The weights of the dual-branch encoder are trained independently, which can learn the feature distribution of different sensor data in a targeted manner. By independently optimizing the weights of the two branches, the problem of accuracy reduction caused by the difference in features of multi-sensor data can be solved. The two modes can be flexibly switched according to the actual data source.

[0032] The beneficial effects of this invention are as follows: 1. Improve the accuracy of mine settlement monitoring: The core adopts KPConv with a rigid core design, which can stably capture the local geometric features of sparse and high-noise point clouds in open-pit mines. Combined with multi-scale downsampling to construct feature pyramids, it can accurately identify small-scale micro-topographic settlement. Compared with traditional 2D or manual feature 3D methods, the monitoring accuracy is significantly improved, meeting the high-precision monitoring needs of mines.

[0033] 2. Adapting to multi-sensor data scenarios: SKPConv-Net is designed with two modes: shared weights and non-shared weights. The non-shared weight mode trains dual-branch weights independently, which can specifically adapt to the feature differences of different sensor data such as satellite and LiDAR, avoid the decrease in accuracy in cross-sensor scenarios, and enhance the compatibility of mine data.

[0034] 3. Reduce the cost of real data annotation: By adopting a strategy of pre-training with simulated data and fine-tuning with real data, after pre-training the model based on the mine simulation dataset, only a small amount of real data is needed for fine-tuning to achieve near-full-scale annotation accuracy, which greatly reduces the workload of high-cost real point cloud annotation in mines.

[0035] 4. Simplify the monitoring process and reduce information loss: The raw 3D point cloud is processed directly throughout the process without the need for rasterization or voxelization preprocessing, thus avoiding the loss of key 3D information; the end-to-end architecture completes the processing from data input to settlement classification without the need for intermediate manual operations, thus simplifying the mine settlement monitoring process. Attached Figure Description

[0036] Figure 1 This is a diagram of the SKPConv-Net neural network architecture in this invention.

[0037] Figure 2 This is a schematic diagram of the process of the present invention. Detailed Implementation

[0038] The specific embodiments of the present invention will now be described in detail with reference to examples and accompanying drawings.

[0039] like Figure 2 The method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference includes the following steps: Step 1: Preprocessing of dual-temporal point cloud data for outdoor to underground scenes: Data preprocessing is the foundation for ensuring effective feature extraction by SKPConv-Net. Addressing issues such as dust interference in outdoor to underground point cloud data, inconsistencies in coordinates due to sensor differences, and uneven point cloud density, the following preprocessing steps are designed: 1.1 Noise Removal: Input dual-temporal original 3D point cloud (Time t1, historical data) and (At time t2, current data) A statistical filtering algorithm is used to remove outlier noise points. Specifically, for each point, the average distance of its K nearest neighbors (K=20) is calculated. If the deviation of the point from the average distance exceeds a preset standard deviation threshold (threshold=1.0), it is identified as a noise point and removed. The formula is as follows:

[0040] in, For point To its The average distance between nearest neighbors, The mean of the average distances between all points. The standard deviation of the mean distance This is the standard deviation threshold.

[0041] 1.2 Coordinate System 1: Since the dual-temporal point clouds may originate from different sensors or different acquisition periods, resulting in coordinate reference differences, the Iterative Closest Point (ICP) algorithm is used to register the dual-temporal point clouds to the same coordinate system. The origin is the mine's pre-set benchmark control point, with the X-axis pointing in the main mining direction and the Z-axis representing the elevation direction. The ICP algorithm solves for the optimal rigid transformation matrix by minimizing the Euclidean distance error between corresponding points in the dual-temporal point clouds. ( For rotation matrix, (where the vector is a translation vector), the formula is as follows:

[0042] in, , Corresponding point pairs in the two-phase point cloud are determined through nearest neighbor search. After registration, Convert to Ensure with The spatial locations are consistent.

[0043] 1.3 Input Feature Construction: If the original point cloud contains supplementary information such as LiDAR intensity and RGB, it is fused with 3D coordinates (X,Y,Z) as the input feature of this network to enhance feature discriminativeness; if it only contains coordinate information, the 3D coordinates are directly used as the input feature, and the input feature dimension is reduced. .

[0044] 1.4 Sub-point cloud segmentation: Considering the large scale of point clouds in open-pit mines, directly inputting them into the network would lead to insufficient computing resources. Therefore, the registered dual-temporal point clouds are segmented... and The system divides the data into vertical cylindrical sub-regions, with the cylinder radius set according to the density of the mine's point cloud: 30m for sparse data and 25m for dense data. This ensures that each sub-point cloud contains sufficient surface detail and potential subsidence areas, while avoiding computational inefficiencies caused by excessively large sub-point clouds. Adjacent sub-point clouds maintain a 30% overlap during segmentation to prevent subsidence areas from being truncated and to ensure the continuity of subsequent monitoring results.

[0045] Step 2: SKPConv-Net neural network design, architecture diagram as shown below Figure 1 As shown, the network is divided into shared weight mode and non-shared weight mode. Figure 1 The diagram shows the neural network architecture in the shared weight mode. Removing the dark blue double-headed dashed arrows representing shared weights reveals the non-shared weight mode. The SKPConv-Net neural network employs a dual-branch encoder and decoder architecture. It extracts point cloud features based on kernel-point convolution KPConv, combining nearest neighbor matching feature difference calculation with skip connections to achieve end-to-end settlement change monitoring. The overall network structure includes an input layer, a dual-branch encoder, a feature difference calculation layer, a decoder, and an output layer. The design of each module is as follows: 2.1 Input Layer Module: Receives the preprocessed dual-temporal sub-point cloud features and coordinates from 1); for each segmented sub-point cloud pair... , (k is the sub-point cloud index), the input feature is and ( , (These are the number of points in the two sub-point clouds, respectively), and their corresponding coordinates are... and The data can be directly input into the dual-branch encoder without the need for rasterization or voxelization preprocessing, thus avoiding information loss.

[0046] 2.2 Dual-branch encoder module: The encoder adopts a dual-branch structure, with the two branches having completely identical architectures, each processing... and Each branch contains 5 layers, each consisting of 2 convolutional blocks. Except for the first layer, the first convolutional block of the remaining layers employs a kernel-based convolutional downsampling design with stride. Multi-scale feature extraction is achieved through recursive downsampling, while maintaining robustness to point clouds of varying densities in the mine. KPConv captures local geometric features by defining kernel distribution and correlation functions within the point cloud neighborhood. Its core formula includes two parts: convolution output calculation and kernel function definition. For point clouds... For any point x in the matrix, the convolution result is the weighted sum of the features of all points in the neighborhood of that point, as shown in the formula:

[0047] in It is a three-dimensional neighborhood centered at x with radius R. Neighboring points Input features, For kernel functions; kernel function g is constructed based on preset kernel points and related functions, and it is applied to neighboring points. offset vector relative to center x Its form is (K is the number of core points, To preset the spatial location of the core point, (for learnable weight matrix), correlation function The distance affected by the core point is positively correlated with R, thus ensuring stable feature extraction in sparse point cloud scenarios.

[0048] To achieve multi-scale feature extraction, the encoder performs recursive downsampling through kernel convolution with stride, with each layer's cell size... With convolution radius according to Recursive definition ( (The initial cell size is set based on the point cloud density). Through this design, the encoder progressively reduces the number of points in the point cloud (downsampling rate of 2) and increases the feature dimension (from...) from layer 1 to layer 5. Up to 512 dimensions), the final output is a multi-scale feature matrix of the dual-temporal sub-point cloud: Branch 1 processing Output and Branch 2 processing Output and Meanwhile, the intermediate features of layers 1-4 of the encoder are passed to the decoder through skip connections to supplement the details lost during downsampling.

[0049] 2.3 Feature Difference Calculation Layer Module: Even in areas without settlement, the number of points in a dual-temporal point cloud can mismatch with their spatial location due to differences in sensor scanning angle and point cloud density, making direct calculation of feature differences prone to errors. Therefore, this layer adopts a nearest neighbor matching-based strategy to calculate the differences in the dual-temporal feature matrices. The specific steps are as follows: right Each point in ,exist Find the point with the closest spatial distance. ,Right now ;extract corresponding feature vector and corresponding feature vector The difference between the two is calculated as a characteristic of change, using the following formula:

[0050] in, The feature difference matrix centrally reflects the settlement-related changes in the two-phase point cloud and serves as the core basis for subsequent settlement classification.

[0051] 2.4 Rotation or Deformation Analysis Module

[0052] To accurately quantify the rotation and non-rigid deformation characteristics of dual-temporal point clouds and provide supplementary geometric basis for settlement classification, a rotation or deformation analysis module has been added after feature difference calculation. This module, based on multi-scale features extracted by KPConv and registered coordinate information, completes the quantitative calculation of rotation angle and deformation intensity in two steps.

[0053] First, the rotation parameters are solved, specifically for the feature difference matrix. The selected candidate regions of change are analyzed based on the intensity of feature differences. Sure, Let be the 90th quantile of the feature differences in the dataset, assuming that the displacement of points within the region conforms to the law of rigid body rotation. For each point within the candidate region... Reuse the kNN matching results from feature difference calculation to obtain its... nearest neighbor in Establish the coordinate correspondence. Solve the rotation matrix R and translation vector t using least squares optimization, with the optimization objective being:

[0054] Where M is the number of points in the candidate region, It is an orthogonal rotation matrix. This is a translation vector. To solve, first calculate the coordinates of the center point of the cloud at the two time points. , After eliminating the translation effect, construct the covariance matrix. ,right Perform SVD decomposition Then the rotation matrix is Rotation matrix Convert to Euler angles ( These correspond to rotation angles around the X, Y, and Z axes, respectively, to quantify the overall rotation trend of the candidate region. The formula is as follows:

[0055] Next, deformation strength calculation is performed. For non-rigid deformation regions, a comprehensive evaluation is conducted using normal vector offset and local strain. The normal vectors at corresponding points at times t1 and t2 are calculated. , Define the normal vector offset Simultaneously, a local neighborhood strain tensor is constructed for each point. Given the set of k nearest neighbors, calculate the covariance matrix of the displacement vector. The principal strain is obtained through eigenvalue decomposition. Define the deformation intensity Ultimately, the rotation angle ( ) and deformation intensity Concatenate into geometric feature vectors , and the feature difference matrix The fusion yields an enhanced feature matrix. This provides more comprehensive feature support for subsequent decoder classification.

[0056] 2.5 Decoder Module: The decoder adopts a 5-layer network structure symmetrical to the encoder. Through a combination of nearest neighbor upsampling, skip connection fusion, and unary convolution, it gradually restores the point cloud resolution while fusing intermediate features from the encoder, ensuring that the output point-scale sinking classification results retain sufficient detail. The specific design is as follows: First, nearest neighbor upsampling is performed on the feature difference matrix. Perform upsampling, setting the sampling rate to 2 for each layer, and use nearest neighbor interpolation to adjust the feature dimensions from the downsampled values. Restore to the original sub-point cloud point count First, it ensures that the final classification result corresponds one-to-one with the original point cloud; second, it performs skip connection fusion. After each layer completes upsampling, the intermediate features of the corresponding layer of the encoder are concatenated with the upsampled features to supplement the local details lost during downsampling. The formula for this process is as follows:

[0057] in Represents the features of the j-th layer of the decoder. Represents the encoder's first Intermediate features of the layer For nearest neighbor upsampling operation, The first step involves feature concatenation; the final step is univariate convolution and nonlinear transformation. For each layer of fused features, univariate convolution is used for dimensionality compression and feature transformation, gradually mapping the feature dimension from 512 dimensions to the settlement category dimension (no settlement, slight settlement, moderate settlement, severe settlement). Simultaneously, the ReLU activation function is introduced to enhance the nonlinear expressive power of the features. The formula is as follows:

[0058] in , These are the weight matrix and bias vector of the j-th unary convolution, respectively.

[0059] 2.6 Output Layer Module: The output layer performs softmax normalization on the transform features of the last layer of the decoder to obtain the sedimentation class probability for each point, as shown in the following formula:

[0060] in, The i-th point belongs to the c-th type of settlement ( No settlement. Slight settlement Moderate settlement. The probability of severe subsidence. Let be the feature value of the i-th point and c-th dimension in the last layer of the decoder.

[0061] To address the problem of class imbalance in mine settlement data, a point-level negative log-likelihood loss function (NLL) is used to optimize the network parameters, as shown in the following formula:

[0062] in, Let i be the actual settlement category label for the i-th point. Here, N represents the predicted probability for the corresponding category, and N is the number of points in the sub-point cloud.

[0063] 2.7 Network Operation Modes: To adapt to multi-sensor data from open-pit to underground mining, SKPConv-Net is designed with two operation modes: Shared weight mode: The dual-branch encoder shares the same network weights and hyperparameters, which is suitable for dual-temporal point clouds with the same sensor and density. This mode has fewer parameters and only needs to train one set of encoder weights, resulting in high training efficiency. It can make full use of the feature consistency of data from the same sensor to improve monitoring accuracy. Non-shared weight mode: The weights of the dual-branch encoder are trained independently, which can learn the feature distribution of different sensor data in a targeted manner. By independently optimizing the weights of the two branches, the problem of accuracy reduction caused by the difference in features of multi-sensor data can be solved. The two modes can be flexibly switched according to the actual data source.

[0064] Example 1: The datasets 3DMtach, 3DLoMatch, and mine dual-temporal point clouds (March 2025 and April 2025) were used for training and testing, including 3D coordinates and optional LiDAR intensity features, in LAS format.

[0065] Step 1: Import the data, process the data, and divide it into training set, validation set, and test set.

[0066] Step 2: Set up the model's operating environment and parameters.

[0067] The model evaluation experiment was conducted on an Ubuntu 20.04 system, using PyTorch 1.10.0, Torch-Points3D, and PyTorch Geometric tools, with GPU acceleration enabled. The hardware consisted of an Intel Xeon Gold CPU, an NVIDIA Tesla V100 GPU, 128GB of DDR4 memory, and a 2TB SSD, meeting the requirements for point cloud computing and storage. The model parameters were configured as follows: 50 iterations, 16 batches, and an optimizer learning rate of 0.001, which decayed by a factor of 0.1 every 25 epochs.

[0068] Step 3: Construct a deep learning-based point cloud change monitoring method. Model structure: By constructing a dual-branch encoder and decoder architecture, end-to-end settlement change monitoring of the original 3D point cloud is achieved. A training strategy of pre-training with simulated data and fine-tuning with real data is adopted to reduce the dependence on real data annotation. At the same time, two network operation modes are used to adapt to multi-sensor data scenarios, ultimately realizing the classification of settlement changes in the dual-temporal point cloud to be monitored.

[0069] The surface subsidence monitoring method for open-pit to underground mining based on point cloud dual-temporal difference proposed in this invention uses deep learning technology to achieve registration between point clouds, outputs surface subsidence data and levels, has high accuracy and robustness, can provide protection for mine safety, and has good application value.

Claims

1. A method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference, characterized in that, The steps are as follows: Step 1: Preprocess the dual-temporal point cloud data of the open-air to underground scene, including denoising, coordinate unification, feature construction, and sub-point cloud segmentation. Step 2: Construct the SKPConv-Net neural network: Adopt a dual-branch encoder and decoder architecture, extract point cloud features based on kernel convolution KPConv, and combine nearest neighbor matching feature difference calculation and jump connection fusion to realize end-to-end settlement change monitoring. The overall network structure includes an input layer, a dual-branch encoder, a feature difference calculation layer, a decoder and an output layer. Step 3: Input the preprocessed dual-temporal point cloud into the trained SKPConv-Net neural network. After feature extraction, difference calculation, deformation analysis and decoding fusion, the displacement of each point and several classification results of no settlement, slight settlement, moderate settlement and severe settlement are obtained, so as to complete the accurate monitoring of surface settlement of open-pit to underground mining.

2. The method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference as described in claim 1, characterized in that, In step 1, the specific method is as follows: Step 1.1 Noise Removal: Input the original 3D point cloud at time t1 (dual-temporal phase). Point cloud at time t2 A statistical filtering algorithm is used to remove outlier noise points. Specifically, for each point, the average distance to its K nearest neighbors is calculated. If the deviation of the point from the average distance exceeds a preset standard deviation threshold, it is identified as a noise point and removed. The formula is as follows: in, For point To its The average distance between nearest neighbors, The mean of the average distances between all points. The standard deviation of the mean distance The standard deviation threshold; Step 1.2 Coordinate System 1: Use the Iterative Closest Point (ICP) algorithm to register the two-phase point clouds to the same coordinate system, with the mine's pre-set benchmark control point as the origin, the X-axis pointing in the main mining direction, and the Z-axis in the elevation direction; solve for the optimal rigid transformation matrix by minimizing the Euclidean distance error between corresponding points in the two-phase point clouds. The formula is as follows: in, Let be a rotation matrix. It is a translation vector. , For corresponding point pairs in the dual-temporal point cloud, the nearest neighbor search is used to determine them. Convert to Ensure with Consistent in spatial location; Step 1.3 Input Feature Construction: If the original point cloud contains supplementary information, it is fused with 3D coordinates (X,Y,Z) as the input features of this network to enhance feature discriminative power; if it only contains coordinate information, the 3D coordinates are directly used as the input features, and the input feature dimension is reduced. ; Step 1.4 Sub-point cloud segmentation: Segment the registered dual-temporal point clouds and The data is divided into vertical cylindrical sub-regions, with the cylinder radius set according to the density of the mine point cloud. The radius is set to 30m for sparse data and 25m for dense data. During the division, adjacent sub-point clouds retain a 30% overlap rate.

3. The method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference as described in claim 1, characterized in that, In step 2, the specific modules of the SKPConv-Net neural network are as follows: Input layer module: Receives the preprocessed dual-temporal sub-point cloud features and coordinates from step 1; for each segmented sub-point cloud pair... , Where k is the sub-point cloud index, and the input feature is and ,in , These represent the number of points in the two sub-point clouds, with corresponding coordinates as follows: and Direct input to the dual-branch encoder without the need for rasterization or voxelization preprocessing; Dual-branch encoder module: The encoder adopts a dual-branch structure, with the two branches having completely identical architectures, each processing... and Each branch contains 5 layers, each consisting of 2 convolutional blocks. Except for the first layer, the first convolutional block of the remaining layers uses a kernel-based convolutional downsampling design with stride, achieving multi-scale feature extraction through recursive downsampling. KPConv captures local geometric features by defining kernel distribution and correlation functions within the point cloud neighborhood. The core formula includes two parts: convolution output calculation and kernel function definition. For point clouds... For any point x in the matrix, the convolution result is the weighted sum of the features of all points in the neighborhood of that point, as shown in the formula: in It is a three-dimensional neighborhood centered at x with radius R. Neighboring points Input features, For kernel functions; kernel function g is constructed based on preset kernel points and related functions, and it is applied to neighboring points. offset vector relative to center x Its form is K is the number of core points. To preset the spatial location of the core point, For the learnable weight matrix, the correlation function in, The influence distance of the core point is positively correlated with R; To achieve multi-scale feature extraction, the encoder performs recursive downsampling through kernel convolution with stride, with each layer's cell size... With convolution radius according to Recursive definition, The initial cell size is set according to the point cloud density. Through this design, the encoder gradually reduces the number of point clouds and increases the feature dimension from layer 1 to layer 5, ultimately outputting a multi-scale feature matrix of dual-temporal sub-point clouds: Branch 1 processing. Output and Branch 2 processing Output and Meanwhile, the intermediate features of layers 1-4 of the encoder are passed to the decoder through skip connections to supplement the details lost during downsampling; Feature difference calculation layer module: This layer adopts a nearest neighbor matching-based strategy to calculate the difference between the two temporal feature matrices. The specific steps are as follows: right Each point in ,exist Find the point with the closest spatial distance. ,Right now ;extract Corresponding feature vector and Corresponding feature vector The difference between the two is calculated as a characteristic of change, using the following formula: in, The feature difference matrix centrally reflects the settlement-related changes in the two-phase point cloud and serves as the core basis for subsequent settlement classification. Rotation or Deformation Analysis Module: A rotation or deformation analysis module has been added after feature difference calculation. This module is based on the multi-scale features extracted by KPConv and the registered coordinate information, and completes the quantitative calculation of rotation angle and deformation intensity in two steps. First, the rotation parameters are solved, specifically for the feature difference matrix. The selected candidate regions of change are analyzed based on the intensity of feature differences. Sure, The 90th quantile of the feature differences in the dataset is given, assuming that the displacement of points within the region conforms to the law of rigid body rotation; for each point within the candidate region... Reuse the kNN matching results from feature difference calculation to obtain its... nearest neighbor in Establish the coordinate correspondence; solve for the rotation matrix R and translation vector t using least squares optimization, with the optimization objective being: Where M is the number of points in the candidate region, It is an orthogonal rotation matrix. This is a translation vector; the solution first calculates the coordinates of the center point of the cloud at two time points. , After eliminating the translation effect, construct the covariance matrix. ,right Perform SVD decomposition Then the rotation matrix is ; rotate matrix Convert to Euler angles ( These correspond to rotation angles around the X, Y, and Z axes, respectively, to quantify the overall rotation trend of the candidate region. The formula is as follows: Next, deformation strength calculation is performed. For non-rigid deformation regions, a comprehensive evaluation is conducted using normal vector offset and local strain. The normal vectors at corresponding points at times t1 and t2 are calculated. , Define the normal vector offset Simultaneously, a local neighborhood strain tensor is constructed for each point. Given the set of k nearest neighbors, calculate the covariance matrix of the displacement vector. The principal strain is obtained through eigenvalue decomposition. Define the deformation intensity Ultimately, the rotation angle ( ) and deformation intensity Concatenate into geometric feature vectors , and the feature difference matrix The fusion yields an enhanced feature matrix. ; Decoder Module: The decoder employs a 5-layer network structure symmetrical to the encoder. Through a combination of nearest neighbor upsampling, skip connection fusion, and unary convolution, it gradually restores the point cloud resolution while fusing intermediate features from the encoder. The first step is nearest neighbor upsampling, targeting the feature difference matrix... Perform upsampling, setting the sampling rate to 2 for each layer, and use nearest neighbor interpolation to adjust the feature dimensions from the downsampled values. Restore to the original sub-point cloud point count Secondly, there's skip connection fusion. After upsampling at each layer, the intermediate features of the corresponding layer from the encoder are concatenated with the upsampled features to supplement the local details lost during downsampling. The formula for this process is: in, Represents the features of the j-th layer of the decoder. Represents the encoder's first Intermediate features of the layer For nearest neighbor upsampling operation, The first step involves feature concatenation; the final step is univariate convolution and nonlinear transformation. For each layer of fused features, univariate convolution is used for dimensionality compression and feature transformation, gradually mapping the feature dimension from 512 dimensions to the settlement category dimension (no settlement, slight settlement, moderate settlement, severe settlement). Simultaneously, the ReLU activation function is introduced to enhance the nonlinear expressive power of the features. The formula is as follows: in , These are the weight matrix and bias vector of the j-th unary convolution, respectively; Output layer module: The output layer performs softmax normalization on the transform features of the last layer of the decoder to obtain the sedimentation class probability of each point, as shown in the following formula: in, Let be the probability that the i-th point belongs to the c-th type of settlement; No settlement. Slight settlement Moderate settlement. Severe settlement Let c be the feature value of the i-th point in the last layer of the decoder; The network parameters are optimized using the point-level negative log-likelihood loss function (NLL), as shown in the following formula: in, Let i be the actual settlement category label for the i-th point. Here, N represents the predicted probability for the corresponding category, and N is the number of points in the sub-point cloud.

4. The method for monitoring surface subsidence in open-pit to underground mining based on point cloud dual-temporal difference as described in claim 1, characterized in that, In step 2, the SKPConv-Net neural network operates in two modes, which can be flexibly switched according to the actual data source: Shared weight mode: The dual-branch encoder shares the same network weights and hyperparameters, which is suitable for dual-temporal point clouds with the same sensor and density. This mode has fewer parameters and only needs to train one set of encoder weights, resulting in high training efficiency. It can make full use of the feature consistency of data from the same sensor to improve monitoring accuracy. Non-shared weight mode: The weights of the dual-branch encoder are trained independently, which can learn the feature distribution of different sensor data in a targeted manner. By independently optimizing the weights of the two branches, the problem of accuracy reduction caused by the difference in features of multi-sensor data can be solved.

Images