A mine area ground surface subsidence monitoring method and system based on unmanned aerial vehicle photogrammetry

Abstract

The application discloses a kind of mine surface subsidence monitoring method and system based on unmanned aerial vehicle photogrammetry, it is related to photogrammetry technical field, method includes: according to unmanned aerial vehicle cruise detection obtains geometric information and radiation information and completes pixel level alignment, constructs joint observation matrix;Rely on subsidence monitor, based on the joint observation matrix of two monitoring nodes carries out curvature change analysis under the decision of decomposition mode and deformation component, obtains curvature flow field vector diagram;According to the vector diagram, mine real-time monitoring and subsidence precursor risk signal prediction are carried out, by calibrating critical streamline and analyzing its time series anomaly pattern, complete subsidence precursor determination and monitoring early warning management.The application solves the technical problems that traditional monitoring means is insufficient in precision for mine surface subsidence monitoring, is difficult to accurately separate interference, cannot identify subsidence precursor in advance, achieves the technical effect of high-precision, global monitoring and active risk warning of surface subsidence under complex open environment.

Description

Technical Field

[0001] This invention relates to the field of photogrammetry technology, and in particular to a method and system for monitoring surface subsidence in mining areas based on UAV photogrammetry. Background Technology

[0002] Surface subsidence monitoring in mining areas is directly related to mine production safety, personnel protection, and geological disaster prevention and control, and is a core aspect of mine safety management. Existing technologies mostly employ single-point measurements and conventional photogrammetry for subsidence detection. While these methods can achieve basic monitoring in stable areas, they are ill-suited to the complex and open spatial characteristics of mining areas. Traditional methods rely solely on single geometric data, failing to integrate multispectral and temperature-related radiation information. Furthermore, they lack computer vision-supported semantic segmentation, deformation feature extraction, and curvature flow field analysis capabilities. This makes it impossible to separate true subsidence deformation from interference signals such as vegetation, temperature, and humidity, hindering the accurate identification of minute deformations and subsidence precursors across the entire area. Consequently, the monitoring data is incomplete and lacks precision, failing to meet the needs of refined assessment and proactive risk management of surface subsidence in mining areas. Summary of the Invention

[0003] This application provides a method and system for monitoring surface subsidence in mining areas based on UAV photogrammetry, which solves the technical problems of insufficient accuracy in monitoring surface subsidence in mining areas, difficulty in accurately separating interference, and inability to identify early signs of subsidence using traditional monitoring methods.

[0004] The first aspect of this application provides a method for monitoring surface subsidence in mining areas based on UAV photogrammetry. The method includes: acquiring geometric and radiometric information based on UAV cruise detection and performing pixel-level alignment to construct a joint observation matrix. The geometric information includes three-dimensional coordinates of a point cloud, and the radiometric information includes band reflectance and temperature inversion values ​​from multispectral images. By developing a subsidence monitor, and based on the joint observation matrix at a determined first and second monitoring time node, performing curvature change analysis under decision-making of decomposed modal components and deformation components to obtain a curvature flow field vector map. Real-time monitoring and management of the mining area are performed based on surface subsidence data from the curvature flow field vector map. Risk signal prediction based on subsidence precursors is also performed based on the curvature flow field vector map for mining area monitoring and early warning management. Specifically, by calibrating critical streamlines in the curvature flow field vector map, subsidence precursors are determined using the time-series anomaly patterns of the critical streamlines under multi-period monitoring.

[0005] A second aspect of this application provides a mining area surface subsidence monitoring system based on UAV photogrammetry. The system includes: a joint observation matrix construction module, used to acquire geometric and radiometric information based on UAV cruise detection and perform pixel-level alignment to construct a joint observation matrix, wherein the geometric information includes three-dimensional coordinates of point clouds, and the radiometric information includes band reflectance and temperature inversion values ​​of multispectral images; a curvature flow field vector map acquisition module, used to obtain a curvature flow field vector map by developing a subsidence monitor and performing curvature change analysis under decision-making of decomposed modal components and deformation components based on the joint observation matrix at a determined first and second monitoring time node; and a mining area monitoring and early warning management module, used to perform real-time monitoring and management of the mining area based on surface subsidence data in the curvature flow field vector map, and to perform risk signal prediction based on subsidence precursors based on the curvature flow field vector map, wherein subsidence precursors are determined by calibrating critical streamlines in the curvature flow field vector map and using the time-series anomaly pattern of the critical streamlines under multi-period monitoring.

[0006] One or more technical solutions provided in this application have at least the following technical effects or advantages:

[0007] This application acquires surface subsidence-related data by collecting geometric and radiometric information from UAV photogrammetry on the surface of the mining area, and then processes it through pixel alignment, bi-branch signal decomposition, and deformation extraction. It calculates information on changes in surface curvature and deformation trends, and makes judgments based on the calibration of critical streamlines and the results of time-series anomaly analysis. This allows for the accurate identification of surface subsidence deformation and subsidence precursors in the mining area, making surface subsidence monitoring in the mining area more accurate and risk warnings more timely and reliable. It achieves the technical effect of high-precision, full-area monitoring and proactive risk warning of surface subsidence in complex and open environments. Attached Figure Description

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

[0009] Figure 1 This is a schematic flowchart of a method for monitoring surface subsidence in mining areas based on unmanned aerial vehicle (UAV) photogrammetry, provided in an embodiment of this application.

[0010] Figure 2 This is a schematic diagram of the structure of a mining area surface subsidence monitoring system based on UAV photogrammetry provided in an embodiment of this application.

[0011] Figure labeling: Joint observation matrix construction module 1, curvature flow field vector map acquisition module 2, mining area monitoring and early warning management module 3. Detailed Implementation

[0012] This application provides a method and system for monitoring surface subsidence in mining areas based on UAV photogrammetry, which solves the technical problems of insufficient accuracy in monitoring surface subsidence in mining areas, difficulty in accurately separating interference, and inability to identify early signs of subsidence using traditional monitoring methods.

[0013] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0014] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this application described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or server that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or modules not explicitly listed or inherent to such processes, methods, products, or devices.

[0015] Example 1, as Figure 1 As shown, a method for monitoring surface subsidence in mining areas based on UAV photogrammetry is provided, wherein the method includes:

[0016] Based on the UAV cruise detection, geometric and radiometric information are acquired and aligned at the pixel level to construct a joint observation matrix. The geometric information includes the three-dimensional coordinates of the point cloud, and the radiometric information includes the band reflectance and temperature inversion values ​​of the multispectral image.

[0017] Specifically, the first step is to complete the flight path planning for the mining area monitoring area. For the mining area surface subsidence monitoring scenario, the cruising altitude is set to 50m to 200m, the forward overlap is no less than 70%, and the lateral overlap is no less than 60%. A UAV equipped with a visible light stereo imaging module, a multispectral imaging module, and a thermal infrared imaging module will perform the cruise detection. With a time synchronization accuracy of no more than 1ms, the imaging acquisition of each module will be triggered synchronously at the preset waypoints to obtain stereo image sequences, multispectral image sequences, and thermal infrared image sequences covering the entire monitoring area.

[0018] Next, using the image control points deployed in the mining area as a reference, aerial triangulation was performed on the acquired stereo image sequence using the bundle method. The image point coordinates of each image were used as the observation values, and collinearity equations were established for each image. Overall adjustment calculations were performed on the exterior orientation elements and densification point coordinates of all images. The adjustment convergence threshold was set to no more than 0.3 pixels. Based on the exterior orientation elements of the adjusted images, a semi-global dense matching algorithm was used to perform dense matching processing to generate high-density three-dimensional point cloud data covering the monitoring area. The point cloud density was no less than 100 points per square meter, and the elevation accuracy was no more than 5 cm. The three-dimensional coordinate information corresponding to each point cloud was extracted as geometric information.

[0019] For the acquired multispectral image sequences, radiometric calibration was performed using the sensor's factory-issued laboratory calibration coefficients. The grayscale values ​​of the images were converted to radiance values. Based on a standard grayscale reference board deployed within the monitoring area, atmospheric correction was performed using an empirical linear method to establish a linear mapping relationship between radiance values ​​and the true surface reflectance. The surface reflectance data for each pixel corresponding to different spectral bands was then calculated. For the acquired thermal infrared image sequences, radiometric calibration was performed using the sensor's factory-issued laboratory calibration coefficients. The grayscale values ​​of the images were converted to radiance values. Combining the measured surface emissivity, average atmospheric temperature, and atmospheric transmittance parameters of the monitoring area, a single-channel algorithm was used to retrieve the surface temperature. First, the on-board brightness temperature was calculated from the radiance values. Then, the true surface temperature was obtained through atmospheric correction and surface emissivity correction. Finally, the surface temperature data corresponding to each pixel was obtained. The band reflectance data and the retrieved temperature values ​​were integrated as radiometric information.

[0020] Subsequently, using the pixel coordinate system of the stereo image sequence as a reference, a mapping relationship between the three-dimensional point cloud coordinates and the image pixel coordinates was established. Corresponding point matching and spatial registration were performed on the multispectral image sequence and the thermal infrared image sequence. For each image, no less than 20 uniformly distributed corresponding control points were selected to solve the homography matrix. Spatial coordinate unification between images was completed through homography matrix transformation. The registered images were resampled using bilinear interpolation. The registration error was controlled within one pixel, so that the three-dimensional coordinates corresponding to the geometric information and the band reflectance and temperature inversion values ​​corresponding to the radiation information were in one-to-one correspondence at the same pixel position, thus completing pixel-level alignment.

[0021] Finally, using pixels with pixel-level alignment as the basic unit, a two-dimensional matrix structure was constructed according to the row and column numbers of the image in the monitoring area. Each row of the two-dimensional matrix corresponds to all pixels with the same row number in the registered image, corresponding to continuous spatial ground units arranged along the flight direction of the UAV within the monitoring area. The total number of rows in the matrix is ​​exactly the same as the total number of rows in the registered image. Each column of the two-dimensional matrix corresponds to all pixels with the same column number in the registered image, corresponding to continuous spatial ground units arranged perpendicular to the flight direction of the UAV within the monitoring area. The total number of columns in the matrix is ​​exactly the same as the total number of columns in the registered image. The geometric and radiometric information of the corresponding pixel is stored synchronously within each matrix unit. Invalid pixels are filled with the neighborhood mean, forming a joint observation matrix covering the entire monitoring area, providing a unified multi-source data foundation for subsequent settlement monitoring and analysis.

[0022] By developing a settlement monitor, and based on the joint observation matrix of the determined first and second monitoring time nodes, a curvature change analysis under the decision of decomposition modal components and deformation components is performed to obtain a curvature flow field vector diagram.

[0023] In this embodiment, the settlement monitor is a software algorithm processing module that integrates signal decomposition, deformation analysis and curvature field calculation functions to process multi-source monitoring data from UAVs and identify surface settlement in mining areas.

[0024] Optionally, a settlement monitor consisting of two parallel data branches with a post-curvature analysis branch is constructed. Using a joint observation matrix of two time series as input, the first branch obtains two types of modal components through adaptive non-recursive signal decomposition, and the second branch extracts deformation components through semantic segmentation and multi-scale morphological analysis. The decomposed modal components and deformation components of the first and second monitoring time nodes are obtained and mapped, respectively. The average curvature, Gaussian curvature and maximum principal curvature of each point are calculated, and corresponding curvature flow field vector diagrams are constructed. The two are then fused to obtain the final curvature flow field vector diagram. The above steps will be explained in detail in the following content.

[0025] Real-time monitoring and management of the mining area are carried out based on the surface subsidence data in the curvature flow field vector diagram. Risk signal prediction based on subsidence precursors is also carried out based on the curvature flow field vector diagram for mining area monitoring and early warning management. Specifically, by marking critical streamlines in the curvature flow field vector diagram, subsidence precursors are determined by the time series anomaly pattern of the critical streamlines under multi-period monitoring.

[0026] In one embodiment of this application, firstly, based on the constructed curvature flow field vector diagram, the core data of surface subsidence contained therein is extracted, including the subsidence location, subsidence rate, and cumulative subsidence of each pixel. A real-time monitoring and management panel for surface subsidence in the mining area is constructed using visualization technology. The subsidence data is synchronously mapped to the mining area topographic map on the monitoring panel according to spatial location. A graded early warning threshold is set. For example, areas with a cumulative subsidence exceeding 50 mm or a subsidence rate exceeding 0.3 mm / day are marked as yellow warnings, and areas with a cumulative subsidence exceeding 100 mm or a subsidence rate exceeding 0.5 mm / day are marked as red warnings. The monitoring data is updated in real time and early warning information is pushed synchronously, realizing real-time monitoring and dynamic management of surface subsidence in the mining area, and providing basic data support for subsequent risk signal prediction.

[0027] Next, the Jacobian matrix is ​​calculated based on the curvature flow field vector diagram to identify four types of critical points: source point, sink point, saddle point, and vortex point. Streamlines are generated by tracing along the direction of the principal eigenvector and a flow field topology skeleton is constructed. The skeleton is then subjected to persistent coherence significance screening to calibrate critical streamlines. Time series anomaly pattern analysis is then carried out on multi-period critical streamlines with a set of preceding and following monitoring time nodes as a single period. Finally, risk signals based on subsidence precursors are generated and used for mine area monitoring and early warning management. This step will be explained in detail in the following sections.

[0028] Furthermore, the method provided in this application embodiment includes:

[0029] Using the time-updated joint observation matrix as input, adaptive non-recursive signal decomposition as a means, and decomposed modal components as output, a first data branch is deployed, wherein the decomposed modal components include physical motion-dominated and biophysical process-dominated components. Using the time-updated joint observation matrix as input, semantic segmentation and multi-scale morphological analysis as means, and deformation components corresponding to the deformation analysis object as output, a second data branch is deployed. The first data branch and the second data branch are paralleled, and a curvature analysis branch is placed afterward as a settlement monitor.

[0030] Specifically, in the process of constructing the settlement monitor, the first data branch is deployed. The first data branch uses the time-updated joint observation matrix as the input data source and adopts adaptive non-recursive signal decomposition as the core processing method to perform feature decomposition processing on the input multi-source time-series observation data. Finally, it outputs the decomposed modal components corresponding to the monitoring area. The decomposed modal components include two categories: components dominated by physical motion and components dominated by biophysical processes, thus completing the classification and decomposition of the temporal change signal of the surface in the mining area.

[0031] Next, the second data branch is deployed. The second data branch uses the joint observation matrix updated in time series, which is from the same source as the first data branch, as the input data source. It uses semantic segmentation and multi-scale morphological analysis as the core processing methods to perform spatial deformation feature extraction processing on the input multi-source observation data. Finally, it outputs the deformation components corresponding to the deformation analysis objects in the monitoring area, thus completing the effective extraction of the spatial deformation features of the surface of the mining area.

[0032] After deploying the two data branches, the first and second data branches are configured to run in parallel. The two data branches synchronously receive the joint observation matrix updated from the same time series source, synchronously execute their respective data processing flows, and synchronously output the corresponding decomposed modal components and deformation components, achieving two-dimensional parallel processing of the same monitoring data.

[0033] Finally, a curvature analysis branch is added at the output of the parallel first and second data branches. This curvature analysis branch takes the decomposed modal components and deformation components output synchronously from the two parallel data branches as input, and carries over the curvature change analysis and processing flow to the subsequent steps, thus completing the overall architecture of the settlement monitor.

[0034] Furthermore, the method provided in this application embodiment includes:

[0035] By performing a Hilbert transform on the original signal in the joint observation matrix, an analytic signal is obtained to construct a single-sided spectrum containing only positive frequency components. The single-sided spectrum of each mode is modulated to the baseband using a mixing exponential term. With the minimum sum of the estimated bandwidths of the single-sided spectra of each mode as a constraint, the alternating direction multiplier method is used for iterative solution to obtain multiple intrinsic mode components with finite bandwidths. For the intrinsic mode components, a first-class mode component based on physical motion and a second-class mode component based on biophysical processes are divided as the decomposed mode components.

[0036] Optionally, the original time-series signal corresponding to each pixel position is first extracted from the joint observation matrix. This signal is real-valued time-series observation data covering multiple monitoring periods, with a sampling interval of 1 to 7 days and a time-series length of no less than 6 complete sampling periods. Discrete Hilbert Transform (DHFT) is performed on each set of original real-valued time-series signals. A convolution operation is performed between the original signals and a preset DHFT kernel to obtain an imaginary component orthogonal to the original signal. The original real-valued signal is combined with the orthogonal imaginary component to construct the corresponding complex-valued analytic signal, completing the analytical processing of the real signal. This eliminates the inherent spectral symmetry redundancy of the real signal, providing a data foundation for the subsequent accurate extraction of the instantaneous frequency and amplitude features of the signal. A Fast Fourier Transform (FFT) is then performed on the complex-valued analytic signal to remove negative frequency components from the transform result, retaining all positive frequency components within the zero-frequency to Nyquist frequency range. This constructs a single-sided spectrum of the complex signal containing only positive frequency components.

[0037] For the constructed single-sided spectrum, the number of intrinsic mode components to be decomposed is preset to 4 to 8. A corresponding mixing exponent term is set for each mode component to be decomposed. The mixing exponent term is composed of the center frequency parameter of each mode and a complex exponential function. By multiplying the mixing exponent term with the corresponding single-sided spectrum of the mode, the high-frequency single-sided spectrum of each mode is shifted down and modulated to the baseband range near zero frequency, thus completing the basebanding of each mode spectrum. Then, the objective function and constraints for optimization are set. The objective function is set to minimize the sum of the estimated bandwidths of the single-sided spectra of all modes to be decomposed. The constraint is the reconstruction constraint that the superposition of all mode components is completely consistent with the original signal. The alternating direction multiplier method is used to perform iterative solution. During the iteration process, the center frequency, bandwidth and component amplitude of each mode are updated alternately. The preset iteration convergence threshold is 1e-6, and the maximum number of iterations is 1000. The iteration process is terminated when the calculation error between two adjacent iterations is less than the preset convergence threshold or the number of iterations reaches the maximum number of iterations. Multiple sets of independent intrinsic mode components with limited bandwidth are output.

[0038] Finally, for all intrinsic modal components obtained by iterative solution, time series characteristic statistics and spectral characteristic analysis are performed respectively. For each intrinsic modal component, the variation law of its time series curve and the distribution characteristics of its spectral energy are extracted. A frequency boundary threshold is set. Intrinsic modal components with spectral energy concentrated in the low frequency band below the frequency boundary threshold, with an absolute value of the linear correlation coefficient of the cumulative change of the time series curve not less than 0.7, and exhibiting monotonous or piecewise linear continuous cumulative change characteristics are classified as Class I modal components based on physical motion. Intrinsic modal components with spectral energy concentrated in the high frequency band above the frequency boundary threshold, with an absolute value of the correlation coefficient with the periodic change parameters of the natural environment not less than 0.6, and exhibiting periodic fluctuation characteristics are classified as Class II modal components based on biophysical processes. The two types of modal components are combined as decomposed modal components.

[0039] By employing a complete process of analytical processing of discrete Hilbert transform signals with defined parameters, iterative solution using baseband spectral modulation and constrained alternating direction multiplier method, and modal classification based on quantized time-frequency characteristics, the effective separation of physical motion signals related to land subsidence and natural environmental interference signals from the original monitoring signals was achieved, providing reliable feature data support for subsequent accurate analysis of land subsidence.

[0040] Furthermore, the method provided in this application embodiment includes:

[0041] The first type of modal component exhibits a low-frequency, high-energy continuous trend signal in terms of time-frequency characteristics. Its time-series curve shows monotonous or piecewise linear cumulative change characteristics, and its spectral energy is concentrated in the low-frequency band with a stable energy accumulation over time. The second type of modal component exhibits high-frequency, periodic fluctuation characteristics in terms of time-frequency characteristics. Its time-series curve has a seasonal rhythm, and its spectral energy is concentrated in the frequency band corresponding to the vegetation growth cycle and temperature and humidity change cycle. Its energy distribution shows regular strength changes with the seasons.

[0042] Specifically, for each set of intrinsic mode components output from the aforementioned steps, time-frequency feature extraction processing is performed. For the time-series data of each intrinsic mode component, discrete Fourier transform is used to perform spectral analysis, and the power spectral density curve corresponding to that component is calculated. The frequency bands where the energy proportion exceeds 90% in the power spectral density curve are statistically concentrated, and the dominant frequency value of that component is determined. At the same time, for the time-series data of each intrinsic mode component, linear regression analysis is used to fit the changing trend of the time-series curve, and the Pearson linear correlation coefficient between the cumulative time-series change and the time variable, as well as the fluctuation period and fluctuation amplitude parameters of the time-series curve, are calculated.

[0043] A frequency boundary threshold of 0.01 cycles / day is set, corresponding to a complete signal variation period of 100 days. This serves as a fixed boundary between the low-frequency and high-frequency bands. Intrinsic mode components with a dominant frequency value less than or equal to this boundary threshold and a corresponding variation period greater than or equal to 100 days are classified as low-frequency signals; intrinsic mode components with a dominant frequency value greater than this boundary threshold and a corresponding variation period less than 100 days are classified as high-frequency signals. This process completes the initial screening and classification of all intrinsic mode components.

[0044] For intrinsic mode components falling into the low-frequency signal range, feature matching verification of a class of mode components is performed. First, the spectral energy distribution of this component is verified, confirming that over 90% of its spectral energy is concentrated in the low-frequency band below the frequency boundary threshold, and that the energy distribution does not exhibit significant abrupt changes over monitoring time, maintaining stable energy accumulation characteristics. Then, the time-series curve characteristics of this component are verified, confirming that the absolute value of the Pearson linear correlation coefficient between its cumulative time-series change and the time variable is not less than 0.7, and that the time-series curve exhibits monotonically increasing, monotonically decreasing, or piecewise linear continuous cumulative change characteristics without significant periodic fluctuations. Intrinsic mode components that simultaneously satisfy both of the above characteristics are identified as a class of mode components dominated by physical motion.

[0045] For intrinsic mode components falling within the high-frequency signal range, feature matching verification of type II mode components was performed. First, the spectral energy distribution of this component was verified, confirming that over 90% of its spectral energy is concentrated in the high-frequency band above the frequency boundary threshold, and the variation period corresponding to the dominant frequency is 30 to 365 days, matching the vegetation growth cycle, annual temperature and humidity variation cycle, and seasonal alternation cycle of the mining area. Subsequently, the time-series curve characteristics of this component were verified, confirming that its time-series curve exhibits periodic fluctuations synchronized with seasonal alternation, and the absolute value of the Pearson correlation coefficient with the temperature, humidity, and vegetation cover data monitored in the same period in the region is not less than 0.6, showing a clear seasonal rhythm, with the energy distribution exhibiting regular changes in strength with seasonal alternation. Intrinsic mode components that simultaneously satisfy both of the above characteristics are identified as type II mode components dominated by biophysical processes.

[0046] Through standardized time-frequency feature extraction, clear frequency threshold division, and quantified feature matching verification process, the two types of modal components were accurately distinguished and characterized, ensuring the complete separation of effective physical motion signals related to land subsidence from natural environmental interference signals, and providing a pure feature data foundation for subsequent accurate analysis of land subsidence deformation.

[0047] Furthermore, the method provided in this application embodiment includes:

[0048] By semantically segmenting the original signals in the joint observation matrix, non-surface features are identified and eliminated, while the effective target distribution for deformation analysis is retained. Multi-scale morphological analysis is performed on the effective target distribution to extract surface features at different scales, wherein the surface features include convex hull, depression, and curvature variation features. Based on the surface features, the deformation components are determined.

[0049] Specifically, firstly, multi-source feature parameters corresponding to each pixel are extracted from the joint observation matrix, including the elevation value corresponding to the three-dimensional coordinates of the point cloud, the normalized vegetation index calculated from the multispectral red band and near-infrared band, and the surface temperature inversion value, and an independent three-dimensional feature vector is constructed for each pixel.

[0050] A CART decision tree semantic segmentation model was constructed based on surface monitoring scenarios in mining areas. This model uses pixel elevation values, normalized difference vegetation index (NDVI), and surface temperature inversion values ​​as input features, and outputs four categories of classification labels: vegetation, buildings, bare soil, and rock. The model is set to a maximum depth of 8 and a minimum leaf node sample size of 5, using the Gini coefficient as a measure of node splitting impurity. UAV photogrammetry data from different areas and time periods within the mining area were collected. The samples include feature parameters of each category of pixels and manually labeled ground truth class labels. The ratio of samples for vegetation, buildings, bare soil, and rock is 2:2:3:3, with a total sample size of no less than 5000. These samples are divided into training and validation datasets in an 8:2 ratio. A gradient boosting strategy is used to iteratively train the model. In each round of training, the optimal features and splitting thresholds for node splitting are calculated based on the training dataset. The number of iterations is set to 100. After each round of iteration, the classification accuracy is calculated using the validation dataset. Training stops when the accuracy improvement is less than 0.5% for 10 consecutive rounds, and the pre-trained CART decision tree semantic segmentation model is obtained.

[0051] Next, a pre-trained CART decision tree semantic segmentation model was used. Using pixel feature vectors as input, pixel-by-pixel semantic classification was performed on all pixels in the monitored area. Fixed decision rules were pre-defined before classification: pixels with a normalized vegetation index greater than 0.3 were classified as vegetation; pixels whose elevation value changed more than 2 meters from the average elevation value of their surrounding 3×3 neighborhood and whose surface temperature was more than 5°C higher than the average temperature of the surrounding bare soil were classified as buildings; and the remaining pixels were classified as bare soil or rock. Ultimately, all pixels were classified into four fixed semantic categories: vegetation, buildings, bare soil, and rock. The overall accuracy of the classification results was no less than 90%. Based on the semantic classification results, pixels classified as vegetation or buildings were marked as non-ground features and were removed entirely. Pixels classified as bare soil or rock were marked as effective surface units. The spatial distribution and corresponding feature data of all effective surface units were fully preserved, forming the effective target distribution for deformation analysis.

[0052] Then, based on the elevation values ​​of each effective surface unit in the effective target distribution, a regular grid digital elevation model matrix covering the monitoring area is constructed. The grid resolution of the matrix is ​​completely consistent with the pixel resolution of the joint observation matrix. Three sets of square morphological structural elements of different sizes (3×3, 5×5, and 7×7) are set, with each set of structural elements corresponding to a surface analysis scale. For the digital elevation model matrix at each analysis scale, morphological opening and closing operations are performed respectively. The opening operation eliminates local bulge noise at the corresponding scale and extracts the surface depression features at the corresponding scale. The closing operation fills the local depression noise at the corresponding scale and extracts the surface convex hull features at the corresponding scale. For the elevation matrices at each scale after morphological processing, the second-order difference method is used to calculate the surface curvature change features corresponding to each effective pixel in the matrix. The calculation window size is consistent with the size of the structural elements at the corresponding analysis scale, thus obtaining the curvature change feature data at each scale. All convex hull, concave and curvature variation features extracted at different scales are spatially matched according to pixel row and column numbers, and multi-scale feature fusion is completed by using equal weight superposition to form a complete surface feature set.

[0053] Finally, using the surface feature set corresponding to the effective target distribution at the first monitoring time node as the spatial reference, and employing the same pixel-level registration rules as those used during the construction of the joint observation matrix, spatial registration processing of the surface feature set at the second monitoring time node was performed. This ensured that the surface features at the two time nodes were aligned in the same spatial coordinate system, with the registration error controlled within one pixel. The elevation difference, curvature change difference, and morphological changes of convex hull and concave features at the same spatial location at the two monitoring time nodes were calculated pixel by pixel to obtain the initial values ​​of deformation for all effective surface units. Outlier removal was performed on the initial deformation values. Values ​​whose absolute value exceeded three times the mean square error of all effective deformation values ​​were marked as outliers, and these outliers were removed entirely. For the initial values ​​of deformation after removing outliers, spatial interpolation is performed using the inverse distance weighting method. During the interpolation process, the power parameter is set to 2, the search neighborhood is set to a 15×15 pixel window, and the maximum number of sampling points for single pixel interpolation does not exceed 20. A continuous deformation spatial distribution field covering the entire monitoring area is generated. The deformation values ​​and spatial distribution characteristics corresponding to each pixel in the deformation spatial distribution field are determined as deformation components.

[0054] By employing a complete process of semantic segmentation with multi-feature decision trees to remove non-ground features with clear decision rules, extraction of multi-scale morphological surface features with fixed fusion rules, and calculation of deformation variables through temporal feature matching with outlier control, the interference of non-deformation features on the surface of the mining area was effectively eliminated, and deformation features related to surface subsidence were accurately extracted, providing reliable deformation data support for the subsequent accurate analysis of surface subsidence in the mining area.

[0055] Furthermore, the method provided in this application embodiment includes:

[0056] At the first monitoring time node, the first decomposition mode component and the first deformation component are acquired; at the second monitoring time node, the second decomposition mode component and the second deformation component are acquired; the first decomposition mode component and the second decomposition mode component, and the first deformation component and the second deformation component are mapped, and curvature field analysis and prediction are performed to determine the surface subsidence data.

[0057] In this embodiment of the application, the first monitoring time node and the second monitoring time node are the time points at which UAV photogrammetry is carried out twice in the monitoring of surface subsidence in the mining area and the corresponding joint observation matrix is ​​generated.

[0058] Specifically, firstly, based on the joint observation matrices corresponding to the determined first and second monitoring time nodes, the corresponding data processing procedures are executed on the joint observation matrices for both periods to obtain the component parameters corresponding to the monitoring data for the two periods. For the joint observation matrix of the first monitoring time node, the original signal in the matrix is ​​fully processed according to the discrete Hilbert transform, spectral modulation iterative solution, and modal classification method described in the previous steps to obtain the first decomposed modal component. This component includes the physical motion-dominated mode and the biophysical process-dominated mode of the first monitoring time node. At the same time, the joint observation matrix is ​​processed according to the semantic segmentation, multi-scale morphological analysis, and deformation calculation method described in the previous steps to obtain the first deformation component, which fully preserves the deformation values ​​and spatial distribution characteristics of the effective surface units of the first monitoring time node.

[0059] For the joint observation matrix at the second monitoring time node, the same processing methods and parameter settings as at the first monitoring time node are used to ensure consistency and comparability of the data processing between the two periods. Similarly, through the signal decomposition process described above, the original signal in the joint observation matrix at the second monitoring time node is processed to obtain the second decomposed modal component, which includes the two types of modal components corresponding to that time node. Through the deformation extraction process described above, the joint observation matrix is ​​processed to obtain the second deformation component, acquiring the deformation values ​​and spatial distribution characteristics of the effective surface units at the second monitoring time node, thus completing the synchronous acquisition of component data from both periods.

[0060] Subsequently, the first decomposition mode component and the second decomposition mode component, as well as the first deformation component and the second deformation component, were subjected to precise spatial location mapping processing. Based on the cell row and column number of the joint observation matrix, the nearest neighbor matching method was used to associate the first decomposition mode component and the second decomposition mode component, as well as the first deformation component and the second deformation component, corresponding to the same spatial location (the same row and column number cell) one-to-one, ensuring that the spatial locations of the component data in the two periods are completely aligned and the matching error is controlled within 1 cell, thus forming a time-series component mapping dataset.

[0061] Finally, based on the mapped time-series component dataset, curvature field analysis and prediction processing were performed. First, the dominant physical motion modal components and deformation components related to surface subsidence were extracted from the two-phase components. Using the second-order difference method, with a 3×3 pixel calculation window, the curvature change values ​​corresponding to the two-phase components were calculated pixel by pixel to construct the surface curvature field for the two monitoring time nodes. Then, a linear fitting prediction method, such as the least squares method, was used as a basis to fit the surface curvature change trend based on the two-phase curvature change values, predicting the short-term surface curvature change pattern. Combining the difference in deformation between the two phases, the surface subsidence rate and cumulative subsidence of each pixel were calculated. Abnormal data with an absolute subsidence rate exceeding 0.5 mm / day were removed, ultimately determining the surface subsidence data covering the entire monitoring area, including core parameters such as subsidence location, subsidence rate, and cumulative subsidence.

[0062] Through a complete process of simultaneous acquisition of components at two time nodes, precise spatial mapping, curvature field analysis, and linear prediction, effective correlation of monitoring data from two periods and accurate analysis of surface subsidence characteristics were achieved, providing comprehensive and reliable data support for dynamic monitoring and risk prediction of surface subsidence in mining areas.

[0063] Furthermore, the method provided in this application embodiment includes:

[0064] For the mapping group of the first decomposition mode component and the second decomposition mode component, a first curvature flow field vector diagram based on curvature change is constructed by calculating the curvature characteristics of each point, wherein the curvature characteristics include mean curvature, Gaussian curvature and maximum principal curvature; for the mapping group of the first deformation component and the second deformation component, a second curvature flow field vector diagram based on curvature change is constructed by calculating the curvature characteristics of each point; the first curvature flow field vector diagram and the second curvature flow field vector diagram are fused together as the curvature flow field vector diagram.

[0065] In one embodiment, the mapping group of the first decomposition mode component and the second decomposition mode component is first extracted, and the dominant physical motion mode component related to land subsidence is selected. Based on the elevation data corresponding to this mode component, a time-series elevation matrix covering the monitoring area is constructed, ensuring that the grid resolution of the matrix is ​​consistent with the joint observation matrix. The second-order difference method is used, with a fixed calculation window of 3×3 pixels, to calculate the curvature features of each point pixel by pixel. During the calculation, a quadratic surface is fitted based on the elevation data of the pixel's neighborhood. By solving the eigenvalues ​​of the quadratic surface, the two principal curvatures corresponding to the pixel are obtained. The average curvature is obtained by solving the two principal curvatures of the fitted surface of the elevation data of the pixel's neighborhood and taking their arithmetic mean. The Gaussian curvature is the product of the two principal curvatures, and the maximum principal curvature is the larger of the two principal curvatures. The three core curvature features of each pixel—average curvature, Gaussian curvature, and maximum principal curvature—are fully obtained. Using the magnitude of the curvature feature as the vector amplitude and the direction of the maximum principal curvature as the vector direction, a vector direction determination rule is set: when the maximum principal curvature is positive, the vector direction points to the extension direction of the surface protrusion; when the maximum principal curvature is negative, the vector direction points to the extension direction of the surface depression. A first curvature flow field vector map based on curvature change is constructed, clearly presenting the curvature change state of the surface micro-topography during the monitoring period.

[0066] For the mapping group of the first and second deformation components, the same method and parameter settings as those used in constructing the first curvature flow field vector map are adopted to ensure the comparability of the two flow fields. Based on the deformation data corresponding to the two deformation components, a time-series deformation elevation matrix is ​​constructed. The second-order difference method and a 3×3 calculation window are also used. During calculation, a quadratic surface is fitted based on the elevation data of the pixel's neighborhood. By solving the eigenvalues ​​of the quadratic surface, the two principal curvatures corresponding to each pixel are obtained. The average curvature, Gaussian curvature, and maximum principal curvature of each point are calculated pixel by pixel to clarify the differences in curvature characteristics at each point. Following the same vector construction rules, with the curvature characteristic value as the amplitude and the direction of the maximum principal curvature as the vector direction, and adhering to the rule that a positive maximum principal curvature points towards the direction of surface bulge extension and a negative maximum principal curvature points towards the direction of surface depression extension, a second curvature flow field vector map based on curvature changes is constructed to accurately reflect the curvature change patterns of the micro-topography corresponding to the deformation components.

[0067] Finally, the first and second curvature flow field vector maps were fused using a fixed-weighted fusion method. Taking into account the importance of the two flow field characteristics, for example, the weight of the first curvature flow field vector map was set to 0.4, and the weight of the second curvature flow field vector map was set to 0.6. This weight allocation aligns with the technical logic in surface subsidence monitoring where deformation components directly reflect the subsidence state and decomposed modal components assist in characterizing the subsidence trend. For pixels at the same spatial location in both flow fields, the vector amplitudes corresponding to their mean curvature, Gaussian curvature, and maximum principal curvature were weighted and summed. The vector direction was selected based on the flow field with the higher weight. The fusion calculation was performed pixel by pixel. In the fused curvature flow field vector map, the vector direction directly corresponds to the evolution direction of the micro-topography of the mining area, and the vector length corresponds to the evolution rate of the micro-topography, fully presenting the dynamic evolution characteristics of the surface micro-morphology within the monitoring area.

[0068] Through a complete process of constructing two flow fields separately and then fusing them with fixed weights, the curvature characteristics corresponding to the decomposed modal components and deformation components were effectively integrated, accurately characterizing the evolution direction and rate of the micro-topography of the mining area, and providing more comprehensive curvature characteristic support for the in-depth analysis and risk prediction of surface subsidence.

[0069] Furthermore, the method provided in this application embodiment includes:

[0070] For the curvature flow field vector diagram, critical points are identified by calculating the Jacobian matrix. These critical points include source points, sink points, saddle points, and vortex points, representing topological singularities where the flow field structure undergoes abrupt changes. Based on these critical points, streamline tracing is performed along the principal eigenvector direction to generate streamlines connecting different critical points, forming a flow field topological skeleton. The flow field topological skeleton is then subjected to saliency screening based on persistent cohomology to calibrate critical streamlines. By performing time-series-based anomaly pattern analysis on critical streamlines from multiple monitoring periods, risk signals based on settlement precursors are generated, where a set of first and second monitoring time nodes is used as a single-period monitoring.

[0071] Optionally, for the curvature flow field vector diagram, the Jacobian matrix corresponding to each pixel in the flow field is calculated using the finite difference method. The Jacobian matrix is ​​constructed as a 2×2 matrix using the partial derivatives of the curvature vector of that pixel in the x and y directions as elements. By solving the eigenvalues ​​and eigenvectors of the Jacobian matrix, critical points in the flow field are identified. Critical point determination rules are set: when both eigenvalues ​​of the Jacobian matrix are positive, the point is a source point, indicating that the flow field vector diffuses outward from this point; when both eigenvalues ​​are negative, the point is a sink point, indicating that the flow field vector converges towards this point; when one eigenvalue is positive and the other is negative, the point is a saddle point, indicating that the flow field vector undergoes a sudden change in direction at this point; when the eigenvalue is a conjugate complex number with a real part of zero, the point is a vortex point, indicating that the flow field vector is in a rotating state at this point. This method completely identifies all source points, sink points, saddle points, and vortex points in the flow field. These critical points are all topological singularities where the flow field structure undergoes a sudden change, and are the core characterization points of precursors to surface subsidence.

[0072] Next, based on all the critical points identified in the above steps, the principal eigenvector direction corresponding to each critical point is extracted. The Euler method is used to perform streamline tracing along the principal eigenvector direction, with a tracing step size of one pixel. Starting from each critical point, the direction and amplitude of the flow field vector are calculated pixel by pixel. Adjacent pixels are connected sequentially along the vector direction to form continuous streamlines, ensuring that the streamlines always conform to the distribution pattern of the flow field vectors and that there are no abrupt changes in direction. The streamlines connecting different critical points are integrated, and broken or overlapping invalid streamlines are removed to form a flow field topology skeleton that can completely reflect the flow field topology structure, clearly presenting the flow pattern and structural distribution of vectors in the flow field.

[0073] Then, a saliency screening based on persistent coherence is performed on the constructed flow field topology skeleton. First, the Vietoris-Rips complex method is used to extract topological features from the streamlines in the topology skeleton. The specific parameters of the Vietoris-Rips complex are: the neighborhood radius is set to 1 pixel unit, the sampling step size is consistent with the grid step size of the joint observation matrix, and the persistence threshold is set to 0.8. Stable topological structures are screened out and temporary streamlines with noise interference are eliminated to ensure that the fused flow field can accurately reflect the evolution law of the surface micromorphology. Next, each streamline is treated as a topological unit, and its persistence is calculated. The persistence calculation focuses on the duration of the streamline's presence in the flow field and its robustness to noise. Specifically, the persistence value is calculated by measuring the existence interval of the topological features corresponding to the streamline at different scales. The length of the interval is the persistence value of the streamline. A higher persistence value indicates that the streamline exists in the flow field for a longer time, has stronger resistance to noise, and has stronger statistical significance, occupying a stable and dominant position in the flow field structure. A lower persistence value indicates that the streamline is more susceptible to noise interference and belongs to temporary structures or noise. A persistence screening threshold of 0.8 is set. Streamlines with persistence values ​​below this threshold are regarded as noise or temporary structures and are removed. Streamlines with persistence values ​​above the threshold and with clear physical significance are retained and labeled as critical streamlines. These critical streamlines can accurately characterize the core propagation path of precursors to surface subsidence.

[0074] Finally, time-series anomaly pattern analysis is conducted on the critical streamlines monitored in multiple periods. When a single critical streamline gradually extends in space and tends to penetrate the zone, or when multiple critical streamlines show regional anomaly convergence or divergence patterns, it is determined whether the anomaly can be resolved by seasonal or periodic factors. If it cannot be resolved, a risk signal of subsidence precursor is generated. This step will be explained in detail in the following content.

[0075] Through a complete process of real-time monitoring and management of curvature flow field data, identification of Jacobi matrix critical points, construction of topological skeleton through streamline tracing, and screening of persistent cohomology significance, the system has achieved accurate identification and prediction of early signs of surface subsidence risk in mining areas, providing a scientific and reliable basis for decision-making in the prevention and control of surface subsidence risk in mining areas.

[0076] Furthermore, the method provided in this application embodiment includes:

[0077] In a time series, when a single critical streamline gradually extends in space and tends to penetrate the zone, a first judgment result is determined; when multiple critical streamlines exhibit regional abnormal convergence or divergence patterns, a second judgment result is determined; if at least one of the first judgment result and the second judgment result is satisfied, it is determined whether the issue can be resolved through seasonal or periodic factors; if not, the risk signal is generated.

[0078] In one embodiment, multiple periods of mining area monitoring data are first acquired. A set of first and second monitoring time nodes are used as a single period of monitoring. Following the methods described above, curvature flow field construction, critical point identification, streamline tracing, and persistent coherence screening operations are performed on each period of monitoring data to obtain the critical streamline corresponding to each period. This ensures that the spatial coordinate system and calculation parameters of the critical streamlines in multiple periods are completely consistent, forming a critical streamline time series dataset covering multiple monitoring periods, providing time series data support with a unified benchmark for subsequent anomaly pattern analysis.

[0079] Next, spatial variation analysis was performed on individual critical streamlines in the critical streamline time series dataset. Using the Euclidean distance measurement method, the length, extension direction, and spatial coordinates of each critical streamline were calculated period by period, and the spatial morphological differences of the same critical streamline in adjacent periods were compared. A streamline extension judgment threshold was set. When the length of a single critical streamline increased by more than 20% compared to the previous period, and the extension direction gradually tended towards known geological weak zones and goaf boundaries in the mining area, the distance between the streamline endpoint and the penetration zone decreased period by period using spatial overlay analysis. This result was determined as the first judgment result, indicating a gradual development trend of surface subsidence precursors. Specifically, combining geological survey data of the mining area and the distribution of goaf areas, the penetration zone was identified as the main subsidence influence zone around the goaf area. The streamline extension direction must be consistent with the direction of the penetration zone to ensure that the streamline analysis is consistent with the actual geological conditions.

[0080] Regional anomaly pattern analysis was performed on multiple streamlines in a multi-phase critical streamline system. Neighborhood analysis was used to calculate the average spacing between adjacent critical streamlines in each phase. Simultaneously, the curvature value of each streamline was calculated using the second-order difference method, and the changes in spacing and curvature between adjacent phases were compared. The aforementioned curvature change analysis is a small-scale superposition analysis. When the spacing between adjacent streamlines shrinks sharply compared to the previous phase (a shrinkage exceeding 30%), or when the curvature value of a streamline suddenly increases (an increase exceeding 50%), it indicates that the region is experiencing concentrated stress. If this anomalous change cannot be explained by other reasonable factors, it suggests that the region may experience step-like subsidence. When multiple critical streamlines exhibiting the aforementioned spacing contraction or curvature abrupt change appear within the monitoring area, and show a regional anomalous convergence or divergence distribution, this is considered the second judgment result.

[0081] Based on the determined first and second judgment results, a seasonal or periodic factor resolution judgment is performed. Pearson correlation coefficient analysis is used to extract seasonal temperature and humidity variation data, vegetation growth cycle data, and precipitation cycle data for the mining area during the same period, calculating the correlation coefficient between the anomaly pattern and the aforementioned periodic factors. A correlation coefficient judgment threshold is set. When the absolute value of the correlation coefficient is not less than 0.6, it indicates that the anomaly pattern can be resolved by seasonal or periodic factors, and no risk signal needs to be generated. When the absolute value of the correlation coefficient is less than 0.6, and there are no other reasonable factors to explain the anomaly, it is judged that it cannot be resolved by seasonal or periodic factors. At this time, a risk signal based on subsidence precursors is generated, containing core information such as the spatial location of the anomaly area, the type of anomaly pattern, and the magnitude of anomaly changes.

[0082] Finally, based on the generated risk signals and the time-series anomaly patterns of critical streamlines in multiple periods, mine area monitoring and early warning management are implemented. In each period's curvature flow field vector diagram, critical streamlines and anomaly areas are marked with special identifiers. Abnormal streamlines are clearly marked with red solid lines, normal and valid streamlines with blue solid lines, and anomaly areas with yellow blocks. The specific location and severity of the abnormal streamlines are simultaneously marked on the monitoring panel for easy identification and handling by staff. For each critical streamline, a corresponding time-series change curve is added to visually present the streamline's changing trend and clearly show the spatial distribution and temporal evolution of the anomaly pattern. A mine area monitoring and early warning management panel is constructed, classifying risk signals into three levels of warning according to the severity of the anomaly. Warning information and anomaly analysis reports are simultaneously pushed to guide staff in conducting targeted on-site investigations and intensified monitoring, achieving accurate early warning and dynamic management based on subsidence precursors.

[0083] By constructing a multi-phase critical streamline time series dataset, analyzing anomaly patterns of single and multiple streamlines, determining the elimination of periodic factors, and implementing a graded early warning management process, the system has achieved accurate identification and scientific early warning of precursory risk signals for surface subsidence in mining areas, providing reliable technical support for the early prevention and timely handling of surface subsidence risks in mining areas.

[0084] In summary, the method for monitoring surface subsidence in mining areas based on UAV photogrammetry provided in this application has the following technical effects:

[0085] This application acquires multi-period surface monitoring data of the mining area through UAV photogrammetry, constructs a curvature flow field vector diagram, and obtains critical streamlines through processing such as Jacobi matrix calculation, streamline tracing, and persistent coherence screening. By analyzing time-series anomaly patterns and combining them with periodic factors, a subsidence precursor risk signal is generated, thereby accurately identifying potential surface subsidence hazards in the mining area. This makes the surface subsidence monitoring and early warning results in the mining area more accurate and reliable, achieving the technical effect of high-precision, full-area monitoring and proactive risk warning of surface subsidence in complex open environments.

[0086] Example 2, as Figure 2 As shown, based on the same inventive concept as in Embodiment 1 above, this application provides a mining area surface subsidence monitoring system based on UAV photogrammetry, the system comprising:

[0087] The joint observation matrix construction module 1 is used to acquire geometric information and radiation information based on UAV cruise detection and perform pixel-level alignment to construct a joint observation matrix. The geometric information includes the three-dimensional coordinates of the point cloud, and the radiation information includes the band reflectance and temperature inversion values ​​of the multispectral image.

[0088] Curvature flow field vector map acquisition module 2 is used to obtain the curvature flow field vector map by developing a settlement monitor and performing curvature change analysis under the decision of decomposing modal components and deformation components based on the joint observation matrix of the determined first monitoring time node and the second monitoring time node.

[0089] The mining area monitoring and early warning management module 3 is used to perform real-time monitoring and management of the mining area based on the surface subsidence data in the curvature flow field vector diagram, and to perform risk signal prediction based on subsidence precursors based on the curvature flow field vector diagram, thereby performing mining area monitoring and early warning management. Specifically, by marking critical streamlines in the curvature flow field vector diagram, subsidence precursors are determined based on the time series anomaly pattern of the critical streamlines under multi-period monitoring.

[0090] Furthermore, the curvature flow field vector diagram acquisition module 2 is used to perform the following steps:

[0091] Using the time-updated joint observation matrix as input, adaptive non-recursive signal decomposition as a means, and decomposed modal components as output, a first data branch is deployed, wherein the decomposed modal components include physical motion-dominated and biophysical process-dominated components. Using the time-updated joint observation matrix as input, semantic segmentation and multi-scale morphological analysis as means, and deformation components corresponding to the deformation analysis object as output, a second data branch is deployed. The first data branch and the second data branch are paralleled, and a curvature analysis branch is placed afterward as a settlement monitor.

[0092] Furthermore, the curvature flow field vector diagram acquisition module 2 is used to perform the following steps:

[0093] By performing a Hilbert transform on the original signal in the joint observation matrix, an analytic signal is obtained to construct a single-sided spectrum containing only positive frequency components. The single-sided spectrum of each mode is modulated to the baseband using a mixing exponential term. With the minimum sum of the estimated bandwidths of the single-sided spectra of each mode as a constraint, the alternating direction multiplier method is used for iterative solution to obtain multiple intrinsic mode components with finite bandwidths. For the intrinsic mode components, a first-class mode component based on physical motion and a second-class mode component based on biophysical processes are divided as the decomposed mode components.

[0094] Furthermore, the curvature flow field vector diagram acquisition module 2 is used to perform the following steps:

[0095] The first type of modal component exhibits a low-frequency, high-energy continuous trend signal in terms of time-frequency characteristics. Its time-series curve shows monotonous or piecewise linear cumulative change characteristics, and its spectral energy is concentrated in the low-frequency band with a stable energy accumulation over time. The second type of modal component exhibits high-frequency, periodic fluctuation characteristics in terms of time-frequency characteristics. Its time-series curve has a seasonal rhythm, and its spectral energy is concentrated in the frequency band corresponding to the vegetation growth cycle and temperature and humidity change cycle. Its energy distribution shows regular strength changes with the seasons.

[0096] Furthermore, the curvature flow field vector diagram acquisition module 2 is used to perform the following steps:

[0097] By semantically segmenting the original signals in the joint observation matrix, non-surface features are identified and eliminated, while the effective target distribution for deformation analysis is retained. Multi-scale morphological analysis is performed on the effective target distribution to extract surface features at different scales, wherein the surface features include convex hull, depression, and curvature variation features. Based on the surface features, the deformation components are determined.

[0098] Furthermore, the curvature flow field vector diagram acquisition module 2 is used to perform the following steps:

[0099] At the first monitoring time node, the first decomposition mode component and the first deformation component are acquired; at the second monitoring time node, the second decomposition mode component and the second deformation component are acquired; the first decomposition mode component and the second decomposition mode component, and the first deformation component and the second deformation component are mapped, and curvature field analysis and prediction are performed to determine the surface subsidence data.

[0100] Furthermore, the curvature flow field vector diagram acquisition module 2 is used to perform the following steps:

[0101] For the mapping group of the first decomposition mode component and the second decomposition mode component, a first curvature flow field vector diagram based on curvature change is constructed by calculating the curvature characteristics of each point, wherein the curvature characteristics include mean curvature, Gaussian curvature and maximum principal curvature; for the mapping group of the first deformation component and the second deformation component, a second curvature flow field vector diagram based on curvature change is constructed by calculating the curvature characteristics of each point; the first curvature flow field vector diagram and the second curvature flow field vector diagram are fused together as the curvature flow field vector diagram.

[0102] Furthermore, the mining area monitoring and early warning management module 3 is used to perform the following steps:

[0103] For the curvature flow field vector diagram, critical points are identified by calculating the Jacobian matrix. These critical points include source points, sink points, saddle points, and vortex points, representing topological singularities where the flow field structure undergoes abrupt changes. Based on these critical points, streamline tracing is performed along the principal eigenvector direction to generate streamlines connecting different critical points, forming a flow field topological skeleton. The flow field topological skeleton is then subjected to saliency screening based on persistent cohomology to calibrate critical streamlines. By performing time-series-based anomaly pattern analysis on critical streamlines from multiple monitoring periods, risk signals based on settlement precursors are generated, where a set of first and second monitoring time nodes is used as a single-period monitoring.

[0104] Furthermore, the mining area monitoring and early warning management module 3 is used to perform the following steps:

[0105] In a time series, when a single critical streamline gradually extends in space and tends to penetrate the zone, a first judgment result is determined; when multiple critical streamlines exhibit regional abnormal convergence or divergence patterns, a second judgment result is determined; if at least one of the first judgment result and the second judgment result is satisfied, it is determined whether the issue can be resolved through seasonal or periodic factors; if not, the risk signal is generated.

[0106] The surface subsidence monitoring system for mining areas based on UAV photogrammetry provided in this embodiment of the invention can execute the surface subsidence monitoring method for mining areas based on UAV photogrammetry provided in any embodiment of the invention, and has the corresponding functional modules and beneficial effects of the execution method.

[0107] Although this application makes various references to certain modules in the system according to the embodiments of this application, any number of different modules can be used and run on user terminals and / or servers. The various units and modules included are only divided according to functional logic, but are not limited to the above division, as long as the corresponding functions can be achieved; in addition, the specific names of each functional unit are only for easy distinction between each other and are not used to limit the scope of protection of this invention.

[0108] The specific embodiments described above do not constitute a limitation on the scope of protection of this application. Those skilled in the art should understand that various modifications, combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the scope of protection of this application. In some cases, the actions or steps described in this application can be performed in a different order than that shown in the embodiments and still achieve the desired results. Furthermore, the processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

Claims

1. A method for monitoring surface subsidence in mining areas based on unmanned aerial vehicle (UAV) photogrammetry, characterized in that, The method includes: Based on the UAV's cruise detection, geometric and radiometric information are acquired and aligned at the pixel level to construct a joint observation matrix. The geometric information includes the three-dimensional coordinates of the point cloud, and the radiometric information includes the band reflectance and temperature inversion values ​​of the multispectral image. By developing a settlement monitor, and based on the joint observation matrix of the determined first and second monitoring time nodes, a curvature change analysis under the decision of decomposed modal components and deformation components is performed to obtain a curvature flow field vector diagram. Real-time monitoring and management of the mining area are carried out based on the surface subsidence data in the curvature flow field vector diagram. Risk signal prediction based on subsidence precursors is carried out based on the curvature flow field vector diagram. Subsidence precursors are determined by marking critical streamlines in the curvature flow field vector diagram and using the time series anomaly pattern of critical streamlines under multi-period monitoring. The development of the settlement monitoring device includes: Using the time-updated joint observation matrix as input, adaptive non-recursive signal decomposition as a means, and decomposed modal components as output, a first data branch is deployed, wherein the decomposed modal components include physical motion-dominated and biophysical process-dominated components. Using the time-updated joint observation matrix as input, semantic segmentation and multi-scale morphological analysis as means, and the deformation components corresponding to the deformation analysis object as output, a second data branch is deployed; The first and second data branches are parallelized, and a curvature analysis branch is added afterward to serve as a settlement monitor. The decision to perform the decomposition of modal components includes: By performing a Hilbert transform on the original signal in the joint observation matrix, an analytic signal is obtained to construct a single-sided spectrum containing only positive frequency components; The single-sided spectrum of each mode is modulated to the baseband by a mixing exponent term. The sum of the estimated bandwidths of the single-sided spectra of each mode is minimized as a constraint. The alternating direction multiplier method is used to iteratively solve the problem and obtain multiple eigenmode components with finite bandwidth. For the intrinsic mode components, a first-class mode component based on physical motion and a second-class mode component based on biophysical processes are divided as the decomposed mode components. The decision to execute the deformation components includes: By performing semantic segmentation on the original signals in the joint observation matrix, non-ground features are identified and eliminated, while the effective target distribution for deformation analysis is preserved. Multi-scale morphological analysis is performed on the effective target distribution to extract surface features at different scales, wherein the surface features include convex hull, depression and curvature variation features; The deformation component is determined based on the surface features. After performing the decision-making process for decomposing the modal components and deformation components, the following steps are included: At the first monitoring time point, the first decomposition mode component and the first deformation component are obtained; At the second monitoring time point, the second decomposition mode component and the second deformation component are obtained; Mapping the first decomposition mode component with the second decomposition mode component and the first deformation component with the second deformation component, performing curvature field analysis and prediction, and determining the surface subsidence data.

2. The method for monitoring surface subsidence in mining areas based on UAV photogrammetry as described in claim 1, characterized in that, The aforementioned modal components exhibit a continuous trend of low frequency and high energy in terms of time-frequency characteristics. The time-series curves show monotonic or piecewise linear cumulative change characteristics, with spectral energy concentrated in the low-frequency band and the energy distribution maintaining stable energy accumulation over time. The two types of modal components exhibit high-frequency, periodic fluctuation characteristics in time-frequency features. The time-series curves have seasonal rhythms, and the spectral energy is concentrated in the frequency bands corresponding to the vegetation growth cycle and temperature and humidity change cycle. The energy distribution shows regular strength changes with the seasonal alternation.

3. The method for monitoring surface subsidence in mining areas based on UAV photogrammetry as described in claim 1, characterized in that, For the mapping group of the first decomposition mode component and the second decomposition mode component, a first curvature flow field vector diagram based on curvature change is constructed by calculating the curvature characteristics of each point, wherein the curvature characteristics include mean curvature, Gaussian curvature and maximum principal curvature; For the mapping group of the first deformation component and the second deformation component, a second curvature flow field vector diagram based on curvature change is constructed by calculating the curvature characteristics of each point; The first curvature flow field vector diagram and the second curvature flow field vector diagram are combined to form the curvature flow field vector diagram.

4. The method for monitoring surface subsidence in mining areas based on UAV photogrammetry as described in claim 1, characterized in that, Based on the curvature flow field vector diagram, risk signal prediction based on settlement precursors is performed, including: For the curvature flow field vector diagram, critical points are identified by calculating the Jacobian matrix. These critical points include source points, sink points, saddle points, and vortex points, which characterize topological singularities where the flow field structure undergoes abrupt changes. Based on the critical point, streamline tracing is performed along the direction of the principal feature vector to generate streamlines connecting different critical points, thus forming the flow field topology skeleton. The critical streamlines are calibrated by performing saliency screening based on persistent cohomology on the flow field topology framework. By performing time-series-based anomaly pattern analysis on critical streamlines from multiple monitoring periods, risk signals based on settlement precursors are generated, with a set of first and second monitoring time nodes serving as a single monitoring period.

5. The method for monitoring surface subsidence in mining areas based on UAV photogrammetry as described in claim 4, characterized in that, Anomaly pattern analysis based on time series of critical streamlines monitored over multiple periods was conducted, including: In a time series, when a single critical streamline gradually extends in space and tends to penetrate the zone, the first judgment result is determined; When multiple critical streamlines exhibit regional abnormal convergence or divergence patterns, a second judgment result is determined. If at least one of the first judgment result and the second judgment result is satisfied, it is determined whether the risk can be mitigated by seasonal or periodic factors. If not, the risk signal is generated.

6. A surface subsidence monitoring system for mining areas based on unmanned aerial vehicle (UAV) photogrammetry, characterized in that, For implementing the method for monitoring surface subsidence in mining areas based on UAV photogrammetry as described in any one of claims 1-5, the system comprises: The joint observation matrix construction module is used to acquire geometric and radiometric information based on UAV cruise detection and perform pixel-level alignment to construct a joint observation matrix. The geometric information includes the three-dimensional coordinates of the point cloud, and the radiometric information includes the band reflectance and temperature inversion values ​​of the multispectral image. The curvature flow field vector map acquisition module is used to obtain the curvature flow field vector map by developing a settlement monitor and performing curvature change analysis under the decision of decomposing modal components and deformation components based on the joint observation matrix of the determined first and second monitoring time nodes. The mining area monitoring and early warning management module is used to perform real-time monitoring and management of the mining area based on the surface subsidence data in the curvature flow field vector diagram, and to predict risk signals based on subsidence precursors based on the curvature flow field vector diagram, thereby performing mining area monitoring and early warning management. Specifically, by marking critical streamlines in the curvature flow field vector diagram, subsidence precursors are determined by the time series anomaly pattern of the critical streamlines under multi-period monitoring.

Images