Goaf three-dimensional subsidence basin reconstruction system based on multi-source remote sensing image fusion

Abstract

The invention relates to the field of remote sensing image processing and geological disaster monitoring, in particular to a goaf three-dimensional subsidence basin reconstruction system based on multi-source remote sensing image fusion, which comprises a heterogeneous data preprocessing module, a differential geometry-based data fusion module, a three-dimensional subsidence basin reconstruction module, a subsidence dynamic monitoring module and a virtual reality interaction module, the system innovatively introduces a differential geometry theory to solve the problem of heterogeneous fusion of a high-resolution optical satellite image, an SAR satellite image and airborne remote sensing Lidar point cloud data, maps multi-source data to a unified feature space through manifold learning, calculates an adaptive fusion weight through curvature analysis, constructs a multi-scale feature correlation matrix based on a geodesic line distance, and achieves the fusion of the high-resolution optical satellite image, the SAR satellite image and the airborne remote sensing Lidar point cloud data. A high-precision fusion image is generated, three-dimensional reconstruction is carried out in combination with subsidence basin geological parameters and a physical constraint method, it is ensured that a reconstruction result conforms to an actual subsidence physical rule, dynamic monitoring and virtual reality interaction of the subsidence process are achieved, and parameter adjustment and model optimization are supported.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of remote sensing image processing and geological disaster monitoring, and particularly relates to a goaf three-dimensional subsidence basin reconstruction system based on multi-source remote sensing image fusion. BACKGROUND

[0002] The goaf surface subsidence in a coal mine is a common geological disaster in the process of coal mining, which poses a serious threat to surface buildings and infrastructure. Accurate monitoring and prediction of the subsidence process of the goaf is of great significance for disaster prevention and mitigation. Traditional goaf subsidence monitoring methods mainly rely on ground leveling and GPS positioning, which have problems such as high cost, low efficiency, and limited coverage.

[0003] With the development of remote sensing technology, satellite remote sensing data has become an important means for monitoring surface subsidence. At present, optical satellite images, SAR (Synthetic Aperture Radar) satellite images, and airborne remote sensing Lidar (Light Detection and Ranging) point cloud data have been widely used in surface subsidence monitoring. However, each single data source has its own limitations: optical satellite images are greatly affected by weather conditions and cannot penetrate clouds; SAR satellite images have all-weather observation capability, but the spatial resolution is relatively low; Lidar point cloud data provides high-precision elevation information, but the cost of acquisition is high and the coverage is limited.

[0004] In the prior art, most multi-source remote sensing data fusion methods use simple weighted averaging or principal component analysis, which cannot effectively process heterogeneous data obtained by different sensors and cannot fully utilize the advantages of each data source. At the same time, traditional subsidence basin reconstruction methods often lack physical constraints, and the reconstruction results do not conform to the actual subsidence physical laws. In addition, existing systems generally lack dynamic updating and interactive functions, making it difficult to meet the needs of real-time monitoring and early warning.

[0005] Therefore, there is an urgent need for a system that can efficiently fuse multi-source remote sensing data and accurately reconstruct the three-dimensional subsidence basin of the goaf to improve the accuracy and efficiency of subsidence monitoring. SUMMARY

[0006] In view of the problems existing in the prior art, the present application provides a goaf three-dimensional subsidence basin reconstruction system based on multi-source remote sensing image fusion, which innovatively introduces differential geometry theory to solve the problem of heterogeneous data fusion, combines physical constraint methods to reconstruct the three-dimensional subsidence basin, and realizes dynamic monitoring and virtual reality interaction.

[0007] The present application provides a goaf three-dimensional subsidence basin reconstruction system based on multi-source remote sensing image fusion, which includes:

[0008] A heterogeneous data preprocessing module for preprocessing high-resolution optical satellite images, SAR satellite images, and airborne remote sensing Lidar point cloud data;

[0009] a differential geometry-based data fusion module, connected with the heterogeneous data preprocessing module, configured to receive the preprocessed multi-source remote sensing data, map the multi-source remote sensing data to a unified feature space based on manifold theory, calculate adaptive fusion weights through curvature analysis, construct a multi-scale feature correlation matrix using geodesic distance, and generate high-precision fused images;

[0010] a three-dimensional subsidence basin reconstruction module, connected with the differential geometry-based data fusion module, configured to receive the high-precision fused images, combine subsidence basin geological parameters, and reconstruct a three-dimensional subsidence basin model using a physical constraint method;

[0011] a subsidence dynamic monitoring module, connected with the three-dimensional subsidence basin reconstruction module, configured to perform time series analysis on the three-dimensional subsidence basin model to realize dynamic monitoring of the subsidence process;

[0012] a virtual reality interaction module, connected with the three-dimensional subsidence basin reconstruction module, configured to generate an interactive virtual reality scene based on the three-dimensional subsidence basin model, and support parameter adjustment and model optimization.

[0013] Preferably, the differential geometry-based data fusion module comprises:

[0014] a manifold learner, configured to regard the multi-source remote sensing data as a point set located in different manifold spaces, construct local geometric structures, and map each data source to a unified manifold space through global mapping optimization;

[0015] a curvature analyzer, connected with the manifold learner, configured to calculate local curvature features of data points in the unified manifold space, construct a structure saliency index based on the curvature features, evaluate information quality of each data source in different regions, and generate adaptive fusion weights;

[0016] a geodesic calculator, connected with the manifold learner, configured to perform geodesic distance calculation on multi-scale features in the unified manifold space, and construct a feature correlation matrix;

[0017] a multi-scale fusioner, connected with the curvature analyzer and the geodesic calculator, configured to perform adaptive fusion on multi-scale features based on the adaptive fusion weights and the feature correlation matrix, and generate the high-precision fused images.

[0018] Preferably, the manifold learner comprises the following specific implementations:

[0019] a feature extraction unit, configured to extract multi-scale texture features, edge features, spectral features, polarization features, scattering features, elevation features, and surface roughness features from the multi-source remote sensing data;

[0020] A local structure modeling unit, connected with the feature extraction unit, is configured to determine K-neighbors for each data point, construct an adjacency graph to represent local relationships, and apply principal component analysis to determine a local coordinate system;

[0021] A global mapping unit, connected with the local structure modeling unit, is configured to optimize global mapping while maintaining local geometry until global structure changes are less than a preset threshold;

[0022] A manifold alignment unit, connected with the global mapping unit, is configured to identify corresponding feature points in different data sources, calculate transformation relationships between manifolds, and map manifolds of different data sources to a unified target manifold space.

[0023] As preferred, the specific implementation of the curvature analyzer includes:

[0024] A local coordinate system establishment unit is configured to establish a tangent space coordinate system in each local region of the unified manifold space and determine a principal direction and a secondary principal direction;

[0025] A curvature tensor calculation unit, connected with the local coordinate system establishment unit, is configured to calculate a first fundamental form and a second fundamental form of a local region and extract a Gaussian curvature and an average curvature;

[0026] A structure saliency measurement unit, connected with the curvature tensor calculation unit, is configured to define structure saliency of a local region based on the curvature value, give high saliency to a high-curvature region, and give low saliency to a low-curvature region;

[0027] A weight generation unit, connected with the structure saliency measurement unit, is configured to combine the structure saliency and data quality evaluation results to construct an initial weight, apply spatial smoothing processing to ensure continuity of weight changes, and perform normalization processing to ensure that all weights sum to 1.

[0028] As preferred, the specific implementation of the geodesic line calculator includes:

[0029] A multi-scale decomposition unit is configured to construct a Gaussian pyramid and a Laplacian pyramid for each data source and extract structure features and texture features of each scale layer;

[0030] A path construction unit, connected with the multi-scale decomposition unit, is configured to connect two points in a feature space to construct an initial path and discretize the path into a finite node sequence;

[0031] A metric calculation unit, connected with the path construction unit, is configured to calculate a local Riemannian metric for each node on the path, taking into account local curvature characteristics of the manifold at the node;

[0032] A path optimization unit, connected with the metric calculation unit, is configured to take the initial path as a search starting point, consider multiple possible forward directions, and select a direction that can minimize the total geodesic distance;

[0033] A distance integration unit, connected with the path optimization unit, is configured to perform Riemann integration along the optimized path to calculate the geodesic distance between the feature points.

[0034] As preferred, the specific implementation of the multi-scale fusioner includes:

[0035] A correlation matrix construction unit is configured to define a distance-to-correlation mapping function based on the geodesic distance, construct a feature correlation matrix, and ensure the symmetry and semi-positivity of the correlation matrix;

[0036] A feature fusion strategy unit, connected with the correlation matrix construction unit, is configured to design cross-scale feature fusion weights based on the correlation matrix, and give higher fusion weights to high-correlation features;

[0037] An inter-scale information transmission unit, connected with the feature fusion strategy unit, is configured to design an information transmission mechanism from a coarse scale to a fine scale, and process areas with inconsistent structures between scales;

[0038] A multi-scale reconstruction unit, connected with the inter-scale information transmission unit, is configured to reconstruct the image from bottom to top, integrate the information of the current layer and the coarser layer at each scale layer, and dynamically adjust the fusion parameters to optimize the visual quality.

[0039] As preferred, the three-dimensional subsidence basin reconstruction module includes:

[0040] A data processing unit is configured to obtain coal mine coordinate information and surface subsidence data corresponding to the coordinate information to form a coal mine surface subsidence data table, preprocess the coal mine surface subsidence data table, and label and frame the target image on the high-precision fused image according to the coal mine coordinate information;

[0041] A parameter inversion unit, connected with the data processing unit, is configured to invert subsidence parameters through coal mining parameters, including constructing a fitting surface through mining area, subsidence depth, and time parameters;

[0042] A basin model calculation unit, connected with the parameter inversion unit, is configured to obtain a surface subsidence basin model through the subsidence parameters and the coal mine coordinate information;

[0043] A three-dimensional model output unit, connected with the basin model calculation unit, is configured to perform three-dimensional modeling on the surface subsidence basin model to view the three-dimensional simulation model of each position in the subsidence area.

[0044] Preferably, the surface subsidence basin model calculation logic in the basin model calculation unit comprises:

[0045] a terrain reference value, representing a ground elevation value corresponding to the position of the mined-out area center or the surface feature point after coordinate conversion, obtained from Lidar point cloud data;

[0046] a mining area fitting surface, representing a mining area fitting surface value corresponding to the position of the mined-out area center or the surface feature point after coordinate conversion;

[0047] a mining rate fitting surface, representing a mining rate fitting surface value corresponding to the position of the mined-out area center or the surface feature point after coordinate conversion;

[0048] a subsidence rate fitting surface, representing a subsidence rate fitting surface value corresponding to the position of the mined-out area center or the surface feature point after coordinate conversion.

[0049] Preferably, the subsidence dynamic monitoring module comprises:

[0050] a data updating unit for obtaining pre-update subsidence basin model data and constructing a pre-update coal mine subsidence parameter model;

[0051] a parameter inversion unit connected to the data updating unit, for obtaining current subsidence basin model data and performing coal mine subsidence parameter inversion;

[0052] a model fusion unit connected to the parameter inversion unit, for fusing the pre-update coal mine subsidence parameter model and the coal mine subsidence parameter inversion result to construct post-update subsidence basin model data;

[0053] a change analysis unit connected to the model fusion unit, for comparing and analyzing the subsidence parameters and the basin model before and after the update, identifying key change areas, and evaluating the subsidence development trend.

[0054] Preferably, the virtual reality interaction module comprises:

[0055] a three-dimensional scene generation unit for generating a three-dimensional scene of the coal mining area based on the three-dimensional subsidence basin model, displaying the subsidence basin model and the mining process of the coal mining area;

[0056] a parameterized adjustment unit connected to the three-dimensional scene generation unit, for adjusting the mining boundary, the dip angle, the strike, the tendency, and the dip angle of the basin slope;

[0057] a dynamic change display unit connected to the parameterized adjustment unit, for drawing and displaying the dynamic change process of the subsidence basin, using dynamic texture mapping and texture deformation technology for subsidence basin terrain drawing;

[0058] A three-dimensional scene rendering unit, connected with the dynamic change display unit, is used for setting a three-dimensional scene light type as a diffuse reflection shadow mode, setting a light source as a parallel light source, setting a scene as a dynamic update, and making a model dynamically interact by adjusting a light coefficient.

[0059] The beneficial effects of the present application include:

[0060] 1. By introducing the differential geometry theory to process multi-source remote sensing data fusion, the heterogeneous data fusion problem is effectively solved, the advantages of each data source are fully utilized, the fusion accuracy is improved by 30%-45% compared with the traditional method, and the edge preservation and detail retention ability are significantly improved.

[0061] 2. Based on manifold learning, curvature analysis and geodesic distance calculation, a complete nonlinear data fusion framework is constructed, which breaks through the limitations of traditional linear fusion methods and is more in line with the physical characteristics of remote sensing data.

[0062] 3. Combined with the geological parameters and physical constraint method of subsidence basin, three-dimensional reconstruction is carried out, so that the reconstruction result is more in line with the actual subsidence physical law, the prediction accuracy of subsidence range is improved by 35%, and reliable basis is provided for disaster warning.

[0063] 4. Dynamic monitoring and virtual reality interaction of subsidence process are realized, parameter adjustment and model optimization are supported, intuitive visual expression of subsidence process is provided, which is helpful for disaster risk assessment and decision support.

[0064] 5. The early identification ability of subsidence is significantly improved, and millimeter-level initial subsidence that is difficult to be found by traditional methods can be identified, the warning time is advanced by 60-90 days, and sufficient time is provided for emergency disposal. BRIEF DESCRIPTION OF DRAWINGS

[0065] Figure 1 It is a schematic diagram of the overall architecture of the system of the present application;

[0066] Figure 2 It is a structural schematic diagram of the data fusion module based on differential geometry;

[0067] Figure 3 It is a schematic diagram of the implementation process of the manifold learner;

[0068] Figure 4 It is a schematic diagram of the working principle of the curvature analyzer;

[0069] Figure 5 It is a schematic diagram of the calculation process of the geodesic calculator;

[0070] Figure 6 It is a schematic diagram of the working process of the multi-scale fusion device;

[0071] Figure 7A structural schematic diagram of the three-dimensional subsidence basin reconstruction module;

[0072] Figure 8 A workflow schematic diagram of the subsidence dynamic monitoring module;

[0073] Figure 9 A functional schematic diagram of the virtual reality interaction module. DETAILED DESCRIPTION

[0074] Reference Figures 1-9 The specific embodiments of the present application will be described in detail below with reference to the accompanying drawings.

[0075] Reference Figure 1 The goaf three-dimensional subsidence basin reconstruction system for multi-source remote sensing image fusion provided by the present application includes a heterogeneous data preprocessing module 1, a differential geometry-based data fusion module 2, a three-dimensional subsidence basin reconstruction module 3, a subsidence dynamic monitoring module 4, and a virtual reality interaction module 5.

[0076] The heterogeneous data preprocessing module 1 is used for preprocessing high-resolution optical satellite images, SAR satellite images, and airborne remote sensing Lidar point cloud data. These preprocessing operations mainly include geometric correction, radiation correction, noise filtering, and other conventional remote sensing data processing methods, which are well known to those skilled in the art and will not be described here.

[0077] The differential geometry-based data fusion module 2 is connected to the heterogeneous data preprocessing module 1 and is used to receive preprocessed multi-source remote sensing data. Based on manifold theory, the multi-source remote sensing data is mapped to a unified feature space, adaptive fusion weights are calculated through curvature analysis, a multi-scale feature correlation matrix is constructed using geodesic distance, and high-precision fused images are generated. This module is the core innovative point of the present application and will be described in detail later.

[0078] The three-dimensional subsidence basin reconstruction module 3 is connected to the differential geometry-based data fusion module 2 and is used to receive high-precision fused images. Combined with subsidence basin geological parameters, a physical constraint method is used to reconstruct a three-dimensional subsidence basin model. This module combines geological parameters with image information to ensure that the reconstruction result conforms to the actual subsidence physical law.

[0079] The subsidence dynamic monitoring module 4 is connected to the three-dimensional subsidence basin reconstruction module 3 and is used for time series analysis of the three-dimensional subsidence basin model to achieve dynamic monitoring of the subsidence process. This module identifies key change areas and evaluates subsidence development trends by comparing subsidence basin models at different time points.

[0080] The virtual reality interaction module 5 is connected to the 3D subsidence basin reconstruction module 3, and is used to generate interactive virtual reality scenes based on the 3D subsidence basin model, supporting parameter adjustment and model optimization. This module provides an intuitive visualization of the subsidence process, which is helpful for disaster risk assessment and decision support.

[0081] Reference Figure 2 The data fusion module 2 based on differential geometry includes a manifold learner 21, a curvature analyzer 22, a geodesic calculator 23, and a multi-scale fusion unit 24.

[0082] The manifold learner 21 is used to treat multi-source remote sensing data as a set of points located in different manifold spaces, construct local geometric structures, and map each data source to a unified manifold space by optimizing the global mapping.

[0083] The curvature analyzer 22 is connected to the manifold learner 21 to calculate the local curvature features of data points in the unified manifold space, construct a structural saliency index based on the curvature features, evaluate the information quality of each data source in different regions, and generate adaptive fusion weights.

[0084] The geodesic calculator 23 is connected to the manifold learner 21 and is used to calculate geodesic distances for multi-scale features in a unified manifold space and construct a feature correlation matrix.

[0085] The multi-scale fusion unit 24 is connected to the curvature analyzer 22 and the geodesic calculator 23 to adaptively fuse multi-scale features based on adaptive fusion weights and feature correlation matrices to generate high-precision fused images.

[0086] Reference Figure 3 The specific implementation of the manifold learner 21 includes a feature extraction unit 211, a local structure modeling unit 212, a global mapping unit 213, and a manifold alignment unit 214.

[0087] The feature extraction unit 211 is used to extract multi-scale texture features, edge features, spectral features, polarization features, scattering features, elevation features, and surface roughness features from multi-source remote sensing data. These feature extraction methods are common in the field of remote sensing image processing and will not be described in detail here.

[0088] The local structure modeling unit 212 is connected to the feature extraction unit 211. It is used to determine the K nearest neighbors for each data point, construct an adjacency graph to represent local relationships, and apply principal component analysis to determine the local coordinate system. In a preferred embodiment of the invention, the value of K is selected between 10 and 25, with the specific value adaptively adjusted according to the data density. When the data distribution is sparse, the K value is smaller (e.g., 10-15); when the data distribution is dense, the K value is larger (e.g., 20-25). Preferably, the adjacency graph is constructed using the KNN (K-Nearest Neighbors) algorithm based on Euclidean distance, and the principal component analysis is implemented using the singular value decomposition method.

[0089] The global mapping unit 213 is connected to the local structure modeling unit 212, and is used to optimize the global mapping while maintaining the local geometric structure until the change in the global structure is less than a preset threshold. In one embodiment of the present invention, the global mapping optimization adopts the following objective function:

[0090] ,

[0091] in: Let be the objective function, representing the weighted sum of squared distances between the mapped data points; The mapped set of data points is a 3D matrix The total number of data points. The dimensions after mapping; The nearest neighbor count is the number of neighboring points considered for each data point. The weight matrix is ​​a 3D matrix, representing points and points The strength of the local relationship between them; Represents the first data point in the mapped data set. The coordinates of a point are a dimensional vector; Represents the Euclidean distance, i.e. ,in For the first The point at the th Coordinates in the dimension.

[0092] Preferably, the weight matrix The Locally Linear Embedding (LLE) algorithm was used to calculate the global structure change threshold, which was set to 0.001. This means that the optimization was considered convergent when the change in the objective function value between two consecutive iterations was less than 0.001. This threshold was an empirical value determined based on extensive experiments, effectively controlling computational complexity while ensuring convergence accuracy.

[0093] Manifold alignment unit 214 is connected to global mapping unit 213 and is used to identify corresponding feature points in different data sources, calculate the transformation relationship between manifolds, and map the manifolds of each data source to a unified target manifold space. In this invention, the correspondence between feature points is identified using the SIFT (Scale-Invariant Feature Transform) algorithm, and the transformation relationship includes a combination of rigid transformations (rotation, translation, scaling) and non-rigid transformations (local deformation). Preferably, the transformation parameters are solved by minimizing the following energy function:

[0094] ,

[0095] in: Let be the energy function, representing the weighted sum of the transformed error and the regularization term; The transformation function maps the points of the source manifold to the target manifold; For the first manifold on the source manifold A corresponding point is a A dimensional vector, where Let be the dimension of the source manifold; For the first on the target manifold A corresponding point is a A dimensional vector, where Let be the dimension of the target manifold; This represents the number of corresponding point pairs; The regularization term controls the smoothness of the transformation and can be represented as the transformation function. The square integral of the second derivative; The regularization coefficient controls the balance between fitting accuracy and smoothness, with a value range of 0.01-0.1 and an optimal value of 0.05.

[0096] Regularization coefficient The choice is based on empirical values. Smaller Values ​​(such as 0.01) emphasize precise matching of corresponding points and are suitable for situations where the number of corresponding points is large and evenly distributed; larger values... Values ​​(such as 0.1) emphasize the smoothness of the transformation and are suitable for situations where the corresponding points are sparse or unevenly distributed. In most goaf monitoring applications, A value of 0.05 can achieve a good balance.

[0097] Reference Figure 4 The specific implementation of the curvature analyzer 22 includes a local coordinate system establishment unit 221, a curvature tensor calculation unit 222, a structural saliency measurement unit 223, and a weight generation unit 224.

[0098] The local coordinate system establishment unit 221 is used to establish a tangent space coordinate system in each local region of the unified manifold space and determine the principal direction and the secondary principal direction. In a preferred embodiment of the present invention, the local region is a region with a resolution of 1.5 to 3 times around the data point. The principal direction and the secondary principal direction are determined by eigenvector analysis, that is, by performing eigenvalue decomposition on the local covariance matrix. The eigenvector with the largest eigenvalue corresponds to the principal direction, and the eigenvector with the second largest eigenvalue corresponds to the secondary principal direction.

[0099] The curvature tensor calculation unit 222 is connected to the local coordinate system establishment unit 221, and is used to calculate the first and second fundamental forms of the local region, extracting Gaussian curvature and mean curvature. In differential geometry, the first fundamental form describes the metric on the manifold, and the second fundamental form describes the degree of curvature of the manifold. For parametric surfaces... The first and second basic forms can be represented as:

[0100] ,

[0101] ,

[0102] in: This is the first basic form, representing the measurement on a curved surface; This is the second basic form, representing the degree of curvature of the surface; , , The coefficients of the first basic form are, where and Representing surfaces Regarding parameters and The partial derivatives, This represents the vector dot product operation; , , The coefficients of the second basic form are, where , , Represent the surface with respect to the parameters and The second-order partial derivative, The unit normal vector is calculated as follows: ,in This represents the vector cross product operation. This represents the magnitude of the cross product result.

[0103] Gaussian curvature and mean curvature It can be calculated using the following formula:

[0104] ,

[0105] ,

[0106] in: Gaussian curvature, representing the intrinsic curvature of the surface at that point, is expressed in units of Gaussian curvature. ; The mean curvature represents the intrinsic curvature of the surface at that point, and its unit is 1000 ppm. Gaussian curvature is the product of the two principal curvatures of a surface at that point, and mean curvature is the arithmetic mean of the two principal curvatures. (Denominator) The determinant of the coefficients of the first fundamental form must be greater than zero, which is always satisfied in regular parametric surfaces.

[0107] The structural saliency measurement unit 223 is connected to the curvature tensor calculation unit 222 and is used to define the structural saliency of local regions based on curvature values, assigning high saliency to high curvature regions and low saliency to low curvature regions. In this invention, the structural saliency index S is defined as:

[0108] ,

[0109] in: The structural saliency index is a dimensionless value that represents the structural importance of an image region. Gaussian curvature; The mean curvature; Let be the gradient of the Gaussian curvature, representing the rate of change of the Gaussian curvature in space, calculated as follows: ; The gradient of the mean curvature represents the rate of change of the mean curvature in space, and is calculated as follows: ; The weighting coefficients are 0.3, 0.3, 0.2, and 0.2, respectively. These weighting values ​​are empirical values ​​that have been verified through a large number of experiments and can achieve good results in most remote sensing scenarios.

[0110] In the above formula Represents absolute value or vector magnitude. For scalars... and , and Represent its absolute value; for vectors and , and Indicates its modulus length, i.e. , .

[0111] The weight generation unit 224 is connected to the structural saliency measurement unit 223. It is used to construct initial weights by combining structural saliency and data quality assessment results, applying spatial smoothing to ensure the continuity of weight changes, and normalization to ensure that all weights sum to 1. Specifically, for spatial location... Data point at, number Weights of each data source The calculation is as follows:

[0112] ,

[0113] in: For the first Data sources in location The weights are numerical values ​​between 0 and 1; For the first Data sources in location The information quality assessment value is a value between 0 and 1. The total number of data sources, in this embodiment These correspond to optical images, SAR images, and Lidar point cloud data, respectively. Indicates all

[0114] The data source is summed.

[0115] Information quality assessment value Factors such as signal-to-noise ratio, sharpness, and detail retention are considered, and the results are calculated using features such as local information entropy and gradient magnitude. In a preferred embodiment of the invention, The calculation method is as follows:

[0116] ,

[0117] in: For the first Data sources in location The signal-to-noise ratio, normalized to interval; The sharpness metric is calculated using local gradients. The degree of detail retention is calculated using local entropy values.

[0118] To ensure the continuity of weight changes, Gaussian filtering is applied for spatial smoothing.

[0119] ,

[0120] in: These are the smoothed weights; For a local window area, it is usually taken as A square region; Let be the Gaussian kernel function, defined as ; The smoothing parameter controls the width of the Gaussian kernel, with a value range of [value range missing]. - The preferred value is .

[0121] Reference Figure 5 The specific implementation of the geodesic calculator 23 includes a multi-scale decomposition unit 231, a path construction unit 232, a metric calculation unit 233, a path optimization unit 234, and a distance integration unit 235.

[0122] The multi-scale decomposition unit 231 is used to construct Gaussian pyramids and Laplacian pyramids for each data source, extracting structural and texture features at each scale level. In a preferred embodiment of the invention, the number of pyramid layers is set to 3-5, and the downsampling ratio between adjacent layers is 2, meaning that each downsampling reduces the image size to 1 / 4 of its original size.

[0123] The path construction unit 232 is connected to the multi-scale decomposition unit 231 and is used to connect two points in the feature space to construct an initial path, discretizing the path into a finite sequence of nodes. The initial path uses straight lines for connection, and the number of discretized nodes depends on the Euclidean distance between the two points. Typically, 1 to 3 nodes are set per unit distance to ensure that the node spacing is less than the effective range of local linearization.

[0124] The metric computation unit 233 is connected to the path construction unit 232 and is used to compute the local Riemannian metric for each node on the path, considering the local curvature properties of the manifold at the node. In Riemannian geometry, the metric tensor... Describing the local geometry on the manifold, for a parameterized surface, the metric tensor can be represented as:

[0125] ,

[0126] in: To measure a tensor, it is a A symmetric positive definite matrix; , , The coefficients of the first basic form are defined in the same way as above.

[0127] In this invention, considering the anisotropy of the feature space, the metric tensor is adjusted as follows:

[0128] ,

[0129] in: The adjusted metric tensor; For the original metric tensor; The feature gradient vector is a 2-dimensional column vector. ,in It is the characteristic function; for The transpose of is a 2-dimensional row vector; This is a weighting coefficient, ranging from 0.5 to 2.0, with a preferred value of 1.0. (Symbol) Representing the outer product of vectors, the result is a The matrix.

[0130] This adjustment increases the metric in directions of dramatic feature changes, thereby increasing the geodesic distance along those directions and helping to preserve feature boundaries. Weighting coefficients The value of is determined based on experimental results, and a larger value indicates a higher probability of ... lower probability Values ​​that enhance boundary preservation may lead to over-segmentation, while smaller values... The value is the opposite.

[0131] The path optimization unit 234 is connected to the metric calculation unit 233 and is used to select the direction that minimizes the total geodesic distance, taking the initial path as the search starting point and considering multiple possible directions of travel. In this invention, path optimization is implemented using the FastMarching Method, which solves the Eikonal equation. ,in Let it be a distance function. Let be the velocity function, and calculate the shortest path from any point to the starting point.

[0132] Distance integration unit 235 is connected to path optimization unit 234 and is used to perform Riemann integration along the optimized path to calculate the geodesic distance between feature points. For path γ(t), t∈[0,1], the geodesic distance d is calculated as follows:

[0133] ,

[0134] Where: d is the geodesic distance, representing the shortest path length between two points on the manifold, with units consistent with the coordinate system; γ(t) is the path parametric equation, t∈[0,1] is the parameter, and γ(0) and γ(1) represent the starting and ending points of the path, respectively; γ(t) is the tangent vector of the path, i.e., the derivative of γ(t) with respect to parameter t; g(γ(t)) is the metric tensor at point γ(t) on the path; the superscript T indicates the transpose of a vector or matrix. Indicates at point At that point, the tangent vector In measurement The length below.

[0135] In the discrete case, the geodesic distance is approximately:

[0136] ,

[0137] in: The distance is a discrete approximation of the geodesic distance; For the first on the path Each node is a 2-dimensional vector; The total number of nodes; For nodes The metric tensor at that location; This represents the difference vector between adjacent nodes; Calculate the squared distance between adjacent nodes in the Riemannian metric.

[0138] Reference Figure 6 The multi-scale fusion unit 24 is specifically implemented by a correlation matrix construction unit 241, a feature fusion strategy unit 242, an inter-scale information transfer unit 243, and a multi-scale reconstruction unit 244.

[0139] The correlation matrix construction unit 241 is used to define a mapping function between distance and correlation based on geodesic distance, construct a characteristic correlation matrix, and ensure the symmetry and positive semidefiniteness of the correlation matrix. In this invention, characteristic correlation... Distance from geodesic lines The relationship is defined as follows:

[0140] ,

[0141] in: Features and characteristics The correlation between them is a value between 0 and 1, with a larger value indicating a stronger correlation. Features and characteristics Geodesic distance between them; It is a natural exponential function; The scaling parameter controls the rate at which the correlation decays with distance, and its value ranges from 0.5 to 2.0, with an optimal value of 1.0.

[0142] This definition ensures the symmetry of the correlation matrix. ) and semidefiniteness. When two characteristics are exactly the same, , When the two features differ greatly, It's very big. Approaching 0. Scale parameter This determines the sensitivity of the correlation to changes in distance; a larger value indicates a lower sensitivity. The value makes the correlation decay more slowly.

[0143] The feature fusion strategy unit 242 is connected to the correlation matrix construction unit 241, and is used to design cross-scale feature fusion weights based on the correlation matrix, assigning higher fusion weights to highly correlated features. Specifically, for scale... and location Features of the place Its fusion weight The calculation is as follows:

[0144] ,

[0145] in: For the first Data sources in location The weights are numerical values ​​between 0 and 1; For the first Data sources in location The information quality assessment value is a value between 0 and 1. The total number of data sources, in this embodiment These correspond to optical images, SAR images, and Lidar point cloud data, respectively. This indicates that the summation is performed across all data sources.

[0146] Characteristic quality assessment value In a preferred embodiment of the present invention, the calculation is as follows:

[0147] ,

[0148] in: For scale In position The local contrast at that point is normalized to the [0,1] interval; The local information entropy is also normalized to the [0,1] interval.

[0149] The inter-scale information transfer unit 243 is connected to the feature fusion strategy unit 242, and is used to design an information transfer mechanism from coarse to fine scales to handle regions with inconsistent structures between scales. In a preferred embodiment of the invention, a Laplacian pyramid fusion strategy is adopted, combining coarse-scale structural information and fine-scale texture details. For each scale... fusion results The calculation is as follows:

[0150] ,

[0151] in: For scale The fusion result; Features at the current scale; For a coarser-scale fusion result; This is an upsampling operation that upsamples the results from a coarser scale to the resolution of the current scale. The weighting coefficient controls the relative contribution of the current scale features and the upsampling results. Its value ranges from 0.3 to 0.7, with an optimal value of 0.5.

[0152] Upsampling operation In this invention, bicubic interpolation is used to maintain image smoothness. For regions with inconsistent structures, Dynamic adjustments give higher weight to scales with clearer structures. The specific adjustment method is as follows:

[0153] ,

[0154] in: For position Adjusted weighting coefficients; and respectively scale and scale In position The structural significance index at the location; The adjustment coefficient has a range of 0.1-0.5, with a preferred value of 0.3.

[0155] Multi-scale reconstruction unit 244 is connected to inter-scale information transfer unit 243 for bottom-up image reconstruction. Each scale layer integrates information from the current layer and coarser layers, dynamically adjusting fusion parameters to optimize visual quality. The final fusion result... The results of the fusion at the finest scale get:

[0156] ,

[0157] In practical applications, post-processing can be performed on the fusion results to optimize visual quality, such as contrast enhancement and edge sharpening. However, these are conventional image processing techniques and will not be discussed in detail here.

[0158] Reference Figure 7 The three-dimensional subsidence basin reconstruction module 3 includes a data processing unit 31, a parameter inversion unit 32, a basin model calculation unit 33, and a three-dimensional model output unit 34.

[0159] The data processing unit 31 is used to acquire coal mine coordinate information and the corresponding surface subsidence to form a coal mine surface subsidence data form. It preprocesses the coal mine surface subsidence data form and marks and selects targets on a high-precision fused image based on the coal mine coordinate information. In one embodiment of the invention, the preprocessing of the data form includes conventional operations such as coordinate system transformation, outlier detection and processing, and data interpolation.

[0160] The parameter inversion unit 32 is connected to the data processing unit 31 and is used to invert subsidence parameters through coal mine mining parameters, including constructing a fitted surface using parameters such as mining area, subsidence depth, and time. The specific implementation process includes:

[0161] First, a fitting curve is obtained based on the mining area. Then, the sample size with different mining area values ​​is selected from the fitting curve, and the fitting surface of the mining area is calculated. Finally, the mining area parameter model and parameters are obtained by fitting using the least squares method.

[0162] Then, based on the settlement depth, a fitting curve is obtained. The sample size with different settlement depth values ​​is selected from the fitting curve, and the settlement depth fitting surface is calculated. The settlement depth parameter model and parameters are obtained by fitting using the least squares method.

[0163] Finally, a fitting curve is obtained based on the time. The sample size with different mining time values ​​is selected from the fitting curve, and the fitting surface of mining time is calculated. The mining time parameter model and parameters are obtained by fitting with the least squares method.

[0164] The basin model calculation unit 33 is connected to the parameter inversion unit 32, and is used to obtain a surface subsidence basin model through subsidence parameters and coal mine coordinate information. In this invention, the calculation logic of the surface subsidence basin model includes:

[0165] The terrain reference value represents the ground elevation value corresponding to the location of the center of the goaf or the surface feature point after coordinate transformation, which is obtained from Lidar point cloud data;

[0166] The fitted surface of the mining area represents the fitted surface value of the mining area corresponding to the location of the center of the goaf or the feature point on the surface after coordinate transformation.

[0167] The mining rate fitting surface represents the mining rate fitting surface value corresponding to the location of the center of the goaf or the surface feature point after coordinate transformation.

[0168] The settlement rate fitting surface represents the settlement rate fitting surface value corresponding to the location of the center of the goaf or the feature point on the surface after coordinate transformation.

[0169] Based on the above parameters, the surface subsidence basin model can be represented as:

[0170] ,

[0171] in: For position The elevation of the ground at the location, in meters; These are topographic reference values, in meters. The surface fitted to the mining area is in square meters; The surface fitted to the mining rate is given in meters per square meter per day (m / (m²·day)) = 1 / day. The surface fitted to the settlement rate is in days, representing the duration of settlement per unit mining area. This reflects the physical process of goaf subsidence: the larger the mining area, the faster the settlement rate, the longer the settlement duration, and ultimately the greater the amount of subsidence.

[0172] The 3D model output unit 34 is connected to the basin model calculation unit 33, and is used to output the 3D model of the surface subsidence basin to view the 3D simulation model of each location within the subsidence area. The 3D modeling adopts TIN (Triangulated Irregular Network) or grid model to generate a visualized 3D subsidence basin model.

[0173] Reference Figure 8 The subsidence dynamic monitoring module 4 includes a data update unit 41, a parameter inversion unit 42, a model fusion unit 43, and a change analysis unit 44.

[0174] Data update unit 41 is used to acquire the pre-update subsidence basin model data and construct the pre-update coal mine subsidence parameter model. The pre-update model data usually comes from historical monitoring results, and the parameter model includes key parameters such as mining area, subsidence depth, and mining time.

[0175] The parameter inversion unit 42 is connected to the data update unit 41 and is used to acquire the current subsidence basin model data and perform coal mine subsidence parameter inversion. The current model data comes from the latest remote sensing data processing results, and the parameter inversion process is similar to the parameter inversion in the 3D subsidence basin reconstruction module.

[0176] Model fusion unit 43 is connected to parameter inversion unit 42 and is used to fuse the pre-update coal mine subsidence parameter model and the coal mine subsidence parameter inversion results to construct the updated subsidence basin model data. The model fusion adopts a time-series filtering method, comprehensively considering the reliability of historical data and current observations. Preferably, Kalman filtering or Bayesian fusion methods are used.

[0177] The change analysis unit 44 is connected to the model fusion unit 43 and is used to compare and analyze the subsidence parameters and basin model before and after the update, identify key change areas, and assess the subsidence development trend. Change analysis mainly examines changes in subsidence amount, subsidence rate, and subsidence range, identifying abnormal areas through change rates and thresholds. In a preferred embodiment of the invention, the subsidence amount change rate threshold is set to 10%, the subsidence rate change rate threshold is set to 15%, and the subsidence range expansion rate threshold is set to 5%. When the change rate of a certain area exceeds the corresponding threshold, the system marks it as an abnormal area for focused monitoring.

[0178] Reference Figure 9The virtual reality interaction module 5 includes a 3D scene generation unit 51, a parameter adjustment unit 52, a dynamic change display unit 53, and a 3D scene rendering unit 54.

[0179] The 3D scene generation unit 51 is used to generate a 3D scene of the coal mining area based on the 3D subsidence basin model, showcasing the coal mining area subsidence basin model and the mining process. The 3D scene includes elements such as surface topography, underground mineral layers, and buildings, and uses realistic texture mapping to enhance the visual effects.

[0180] The parameter adjustment unit 52 is connected to the three-dimensional scene generation unit 51 and is used to adjust the dip angle, strike, dip direction, and dip angle of the mining boundary and basin slope. In one embodiment of the present invention, the adjustment range of the slope dip angle is 10°-45°, the adjustment range of the strike angle is 0°-360°, and the adjustment range of the dip angle is 0°-90°.

[0181] The dynamic change display unit 53 is connected to the parametric adjustment unit 52 and is used to draw and display the dynamic change process of the subsidence basin. Dynamic texture mapping and texture deformation techniques are used to draw the subsidence basin topography. The dynamic change display is based on time-series data, showing the evolution of the subsidence process on a monthly or quarterly basis. Preferably, an interpolation method is used to generate a continuous change process, improving the smoothness of the animation.

[0182] The 3D scene rendering unit 54 is connected to the dynamic change display unit 53, and is used to set the 3D scene lighting type to diffuse shadow mode, set the light source to parallel light source, set the scene to dynamic update, and enable dynamic interaction of the model by adjusting the lighting coefficient. In a preferred embodiment of the present invention, the lighting direction is set to shine from the northwest, the incident angle is 45°, and the shadow intensity is set to 0.6 to enhance the three-dimensionality and realism of the 3D scene.

[0183] This invention addresses the challenge of heterogeneous data fusion by introducing differential geometry theory and combining it with physical constraint methods to reconstruct three-dimensional subsidence basins, achieving high-precision monitoring and prediction of subsidence in mining areas. In practical applications, the system achieves an accuracy rate of over 85% in identifying initial subsidence at the millimeter level, improves the accuracy of subsidence range prediction by 35%, and advances the warning time by 60–90 days.

[0184] Furthermore, the system's virtual reality interactive function provides an intuitive visualization tool for disaster risk assessment and decision support, helping to improve the efficiency and effectiveness of disaster prevention and mitigation. Through parameter adjustment and model optimization, users can simulate subsidence under different mining schemes, providing a scientific basis for mining planning.

[0185] In summary, the multi-source remote sensing image fusion-based three-dimensional subsidence basin reconstruction system for mining subsidence areas provided by this invention solves the problem of heterogeneous data fusion by innovatively introducing differential geometry theory, thereby improving the accuracy and efficiency of subsidence monitoring. It has significant theoretical value and promising practical application prospects.

[0186] The embodiments described above are merely illustrative of specific implementations of the present invention, and while the descriptions are detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A three-dimensional subsidence basin reconstruction system for mined-out areas based on multi-source remote sensing image fusion, characterized in that, include: The heterogeneous data preprocessing module is used to preprocess high-resolution optical satellite imagery, SAR satellite imagery, and airborne remote sensing LiDAR point cloud data. A data fusion module based on differential geometry is connected to the heterogeneous data preprocessing module. It is used to receive the preprocessed multi-source remote sensing data, map the multi-source remote sensing data to a unified feature space based on manifold theory, calculate adaptive fusion weights through curvature analysis, construct a multi-scale feature correlation matrix using geodesic distance, and generate a high-precision fused image. The three-dimensional subsidence basin reconstruction module is connected to the differential geometry-based data fusion module. It is used to receive the high-precision fused image, combine it with the geological parameters of the subsidence basin, and reconstruct the three-dimensional subsidence basin model using a physical constraint method. The subsidence dynamic monitoring module is connected to the three-dimensional subsidence basin reconstruction module and is used to perform time-series analysis on the three-dimensional subsidence basin model to realize dynamic monitoring of the subsidence process. The virtual reality interaction module is connected to the three-dimensional subsidence basin reconstruction module and is used to generate an interactive virtual reality scene based on the three-dimensional subsidence basin model, supporting parameter adjustment and model optimization.

2. The system according to claim 1, characterized in that, The differential geometry-based data fusion module includes: A manifold learner is used to treat the multi-source remote sensing data as a set of points located in different manifold spaces, construct local geometric structures, and map each data source to a unified manifold space by optimizing the global mapping. A curvature analyzer, connected to the manifold learner, is used to calculate the local curvature features of data points in the unified manifold space, construct a structural saliency index based on the curvature features, evaluate the information quality of each data source in different regions, and generate adaptive fusion weights. A geodesic calculator, connected to the manifold learner, is used to calculate geodesic distances to multi-scale features in the unified manifold space and construct a feature correlation matrix. A multi-scale fusion unit, connected to the curvature analyzer and the geodesic calculator, is used to adaptively fuse multi-scale features based on the adaptive fusion weights and the feature correlation matrix to generate the high-precision fused image.

3. The system according to claim 2, characterized in that, The specific implementation of the manifold learner includes: The feature extraction unit is used to extract multi-scale texture features, edge features, spectral features, polarization features, scattering features, elevation features, and surface roughness features from the multi-source remote sensing data. A local structure modeling unit, connected to the feature extraction unit, is used to determine the K nearest neighbors for each data point, construct an adjacency graph to represent local relationships, and apply principal component analysis to determine the local coordinate system; A global mapping unit, connected to the local structure modeling unit, is used to optimize the global mapping while maintaining the local geometric structure, until the change in the global structure is less than a preset threshold. The manifold alignment unit, connected to the global mapping unit, is used to identify corresponding feature points in different data sources, calculate the transformation relationship between manifolds, and map the manifolds of each data source to a unified target manifold space.

4. The system according to claim 2, characterized in that, The specific implementation of the curvature analyzer includes: A local coordinate system establishment unit is used to establish a tangent space coordinate system in each local region of the unified manifold space and determine the principal direction and secondary principal direction. The curvature tensor calculation unit is connected to the local coordinate system establishment unit and is used to calculate the first and second basic forms of the local region and extract the Gaussian curvature and the average curvature. The structural saliency measurement unit, connected to the curvature tensor calculation unit, is used to define the structural saliency of a local region based on the curvature value, assigning high saliency to high curvature regions and low saliency to low curvature regions. The weight generation unit, connected to the structural saliency measurement unit, is used to construct initial weights by combining the structural saliency and data quality assessment results, apply spatial smoothing to ensure the continuity of weight changes, and ensure that all weights sum to 1 through normalization.

5. The system according to claim 2, characterized in that, The specific implementation of the geodesic calculator includes: Multi-scale decomposition units are used to construct Gaussian and Laplacian pyramids for each data source and extract structural and texture features at each scale level. The path construction unit, connected to the multi-scale decomposition unit, is used to connect two points in the feature space to construct an initial path and discretize the path into a finite sequence of nodes. The metric calculation unit, connected to the path construction unit, is used to calculate the local Riemann metric for each node on the path, taking into account the local curvature characteristics of the manifold at the node. The path optimization unit, connected to the metric calculation unit, is used to take the initial path as the search starting point, consider multiple possible directions of progress, and select the direction that minimizes the total geodesic distance. The distance integration unit, connected to the path optimization unit, is used to perform Riemann integration along the optimized path to calculate the geodesic distance between feature points.

6. The system according to claim 2, characterized in that, The specific implementation of the multi-scale fusion processor includes: The correlation matrix construction unit is used to define a mapping function between distance and correlation based on the geodesic distance, construct a feature correlation matrix, and ensure the symmetry and positive semidefiniteness of the correlation matrix. The feature fusion strategy unit, connected to the correlation matrix construction unit, is used to design cross-scale feature fusion weights based on the correlation matrix and assign higher fusion weights to highly correlated features. The inter-scale information transfer unit, connected to the feature fusion strategy unit, is used to design an information transfer mechanism from coarse scale to fine scale and to handle regions with inconsistent structures between scales. A multi-scale reconstruction unit, connected to the inter-scale information transfer unit, is used to reconstruct images from the bottom up. Each scale layer integrates information from the current layer and coarser layers, and dynamically adjusts the fusion parameters to optimize visual quality.

7. The system according to claim 1, characterized in that, The three-dimensional subsidence basin reconstruction module includes: The data processing unit is used to acquire coal mine coordinate information and the corresponding surface subsidence amount to form a coal mine surface subsidence data form, preprocess the coal mine surface subsidence data form, and mark and select the target image on the high-precision fused image according to the coal mine coordinate information. The parameter inversion unit, connected to the data processing unit, is used to invert subsidence parameters through coal mining parameters, including constructing a fitted surface using parameters such as mining area, subsidence depth and time. The basin model calculation unit is connected to the parameter inversion unit and is used to obtain the surface subsidence basin model through the subsidence parameters and the coal mine coordinate information. The 3D model output unit is connected to the basin model calculation unit and is used to perform 3D modeling and output of the surface subsidence basin model to view the 3D simulation model of each location in the subsidence area.

8. The system according to claim 7, characterized in that, The surface subsidence basin model calculation logic in the basin model calculation unit includes: The terrain reference value represents the ground elevation value corresponding to the location of the center of the goaf or the surface feature point after coordinate transformation, which is obtained from Lidar point cloud data; The fitted surface of the mining area represents the fitted surface value of the mining area corresponding to the location of the center of the goaf or the feature point on the surface after coordinate transformation. The mining rate fitting surface represents the mining rate fitting surface value corresponding to the location of the center of the goaf or the surface feature point after coordinate transformation. The settlement rate fitting surface represents the settlement rate fitting surface value corresponding to the location of the center of the goaf or the feature point on the surface after coordinate transformation.

9. The system according to claim 1, characterized in that, The subsidence dynamic monitoring module includes: The data update unit is used to acquire the data of the subsidence basin model before the update and to construct the coal mine subsidence parameter model before the update. The parameter inversion unit, connected to the data update unit, is used to acquire the current subsidence basin model data and perform coal mine subsidence parameter inversion. The model fusion unit, connected to the parameter inversion unit, is used to fuse the coal mine subsidence parameter model before the update and the coal mine subsidence parameter inversion results to construct the updated subsidence basin model data. The change analysis unit, connected to the model fusion unit, is used to compare and analyze the subsidence parameters and basin model before and after the update, identify key change areas, and assess the subsidence development trend.

10. The system according to claim 1, characterized in that, The virtual reality interaction module includes: The three-dimensional scene generation unit is used to generate a three-dimensional scene of the coal mining area based on the three-dimensional subsidence basin model, and to display the subsidence basin model and mining process of the coal mining area. The parameterized adjustment unit, connected to the three-dimensional scene generation unit, is used to adjust the dip angle, strike, dip direction, and dip angle of the mining boundary and basin slope. The dynamic change display unit, connected to the parameterized adjustment unit, is used to draw and display the dynamic change process of the subsidence basin, and uses dynamic texture mapping and texture deformation technology to draw the subsidence basin terrain. The 3D scene rendering unit, connected to the dynamic change display unit, is used to set the 3D scene lighting type to diffuse shadow mode, set the light source to parallel light source, set the scene to dynamic update, and enable dynamic interaction of the model by adjusting the lighting coefficient.

Images