A method and related equipment for fusing multi-source monitoring data of surface subsidence in mining areas

Abstract

This application provides a method and related equipment for fusing multi-source monitoring data of surface subsidence in mining areas. The method includes: acquiring monitoring image data of a preset study area by an unmanned aerial vehicle (UAV), and generating a digital surface model based on the monitoring image data; cutting the digital surface model into multiple grids, and determining the grid data of each grid based on the digital surface model; deploying monitoring points in at least some of the grids, and acquiring first coordinate data and second coordinate data of the monitoring points; correcting the grid data of each grid based on the first coordinate data and second coordinate data of the monitoring points to obtain corrected grid data; and determining the corrected grid data of all grids as multi-source monitoring fusion data of surface subsidence in mining areas; wherein the accuracy of the second coordinate data is higher than that of the first coordinate data. This data fusion method achieves good spatial and temporal fusion and high accuracy of the fused data.

Description

Technical Field

[0001] This application relates to the field of surface monitoring technology, and in particular to a method and related equipment for fusing multi-source monitoring data of surface subsidence in mining areas. Background Technology

[0002] Monitoring surface subsidence is crucial for geological safety. Traditional single monitoring methods have limitations in data accuracy and application, making it difficult to comprehensively capture dynamic changes. Therefore, multi-source data fusion technology has become an important development direction. However, this technology faces two major challenges: first, data acquired by different monitoring methods vary greatly in dimensionality and spatiotemporal resolution, exhibiting significant geographic spatial heterogeneity, making direct assimilation and unification of the data difficult; second, how to select appropriate methods to fully leverage the advantages of each monitoring method and achieve efficient and accurate multi-source information fusion remains a key issue that urgently needs to be addressed. Summary of the Invention

[0003] In view of this, the purpose of this application is to propose a method and related equipment for fusing multi-source monitoring data of surface subsidence in mining areas, so as to solve the problem that multi-source monitoring data are difficult to assimilate and unify due to differences in dimensionality, spatiotemporal resolution and geographic spatial heterogeneity under the limitations of the accuracy and application of single monitoring methods, and the lack of efficient and accurate fusion methods.

[0004] To achieve the above objectives, this application provides a method for fusing multi-source monitoring data of surface subsidence in mining areas, comprising:

[0005] Acquire monitoring image data of a preset research area by a drone, and generate a digital surface model based on the monitoring image data;

[0006] The digital surface model is cut into multiple grids, and the grid data of each grid is determined based on the digital surface model;

[0007] Monitoring points are deployed in at least a portion of the grid, and first coordinate data and second coordinate data of the monitoring points are acquired;

[0008] Based on the first and second coordinate data of the monitoring point, the grid data of each grid is corrected to obtain the corrected grid data.

[0009] The corrected raster data of all graticules were identified as multi-source monitoring and fusion data of surface subsidence in the mining area;

[0010] The first coordinate data is determined based on the digital surface model, while the second coordinate data is not determined based on the digital surface model, and the accuracy of the second coordinate data is higher than that of the first coordinate data.

[0011] Optionally, the second coordinate data includes coordinate data obtained through manual monitoring and coordinate data obtained through satellite positioning.

[0012] Optionally, the step of correcting the raster data of each grid cell based on the first and second coordinate data of the monitoring point to obtain corrected raster data includes:

[0013] In the pre-constructed improved radial basis function interpolation algorithm model, the raster data of each grid is corrected based on the first coordinate data and the second coordinate data of the monitoring point to obtain the corrected raster data.

[0014] Optionally, the pre-construction process of the improved radial basis function interpolation algorithm model includes the following steps:

[0015] Determine the terrain roughness and tilt of each grid cell;

[0016] Based on the terrain roughness of all the grids, the mean terrain roughness of the study area is determined;

[0017] For each of the aforementioned grid cells, the following steps are performed:

[0018] Determine the shape parameters of a grid based on its terrain roughness and tilt.

[0019] Based on the shape parameters of the grid, the grid data of the grid, and the first coordinate data of the corresponding target monitoring point, a Gaussian function model is constructed; each grid corresponds to at least one target monitoring point, which is a monitoring point located within the grid or a monitoring point located around the grid;

[0020] Based on the grid data of the grid and the first coordinate data of the corresponding target monitoring point, a multi-quadratic function model is constructed;

[0021] Based on the mean terrain roughness and the terrain roughness of the grid, a first weight for the Gaussian function model and a second weight for the multiple quadratic function model are determined.

[0022] The influence factor of the monitoring point on the grid is determined based on the first weight, Gaussian function model, second weight and multiple quadratic function model;

[0023] Determine the elevation residuals of the first and second coordinate data of at least one target monitoring point;

[0024] An improved radial basis function interpolation algorithm model is constructed based on the first coordinate data of at least one target monitoring point, the elevation residual, and the influence factor on the grid.

[0025] Optionally, determining the terrain roughness and tilt of each of the grid cells includes:

[0026] Obtain the raster data of a given grid cell and all its adjacent grid cells. Based on the raster data of the given grid cell and all its adjacent grid cells, obtain the elevation data of the given grid cell and all its adjacent grid cells.

[0027] Based on the elevation data of this grid and all adjacent grids, obtain the average elevation data;

[0028] The terrain roughness of the grid is determined based on the elevation data and average elevation data of the grid and all adjacent grids.

[0029] Optionally, determining the terrain roughness and tilt of each of the grid cells further includes:

[0030] Obtain raster data of a grid cell and its neighboring grid cells in the x-direction. Based on the raster data of the grid cell and its neighboring grid cells in the x-direction, obtain the elevation data of the grid cell and its neighboring grid cells in the x-direction.

[0031] Based on the elevation data of the grid and its adjacent grids in the x-direction, the height gradient of the grid in the x-direction is obtained;

[0032] Obtain raster data of a grid cell and its neighboring grid cells in the y-direction. Based on the raster data of the grid cell and its neighboring grid cells in the y-direction, obtain the elevation data of the grid cell and its neighboring grid cells in the y-direction.

[0033] Based on the elevation data of the grid and its adjacent grids in the y-direction, the height gradient of the grid in the y-direction is obtained;

[0034] Obtain the elevation scaling factor of the digital surface model;

[0035] The tilt of the grid is determined based on the height gradient of the grid in the x-direction and the height gradient in the y-direction, as well as the elevation scaling factor.

[0036] The x-direction and y-direction are located on the same horizontal plane.

[0037] Optionally, determining the mean terrain roughness of the study area based on the terrain roughness of all the grids includes:

[0038] Based on the terrain roughness of all the grids and the mean terrain roughness of the study area, the standard deviation of the terrain roughness of the study area is determined;

[0039] The method of determining the shape parameters of a grid based on its terrain roughness and tilt also includes:

[0040] The shape parameters of the grid are determined based on the mean and standard deviation of the terrain roughness of the study area, as well as the terrain roughness and tilt of the grid.

[0041] Optionally, the influence factor of the monitoring point on the grid is determined based on the first weight, Gaussian function model, second weight, and multiple quadratic function model, expressed by the formula:

[0042] ;

[0043] in, The influence factor of the i-th target monitoring point on the grid; It is a Gaussian function model; It is a model of multiple quadratic functions; The Euclidean distance between the grid and the corresponding target monitoring point on the horizontal plane is calculated using the grid data of the grid and the first coordinate data of the corresponding target monitoring point. It is the first weight; is the second weight, and c is the rate constant.

[0044] Based on the same inventive concept, this disclosure also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable by the processor, wherein the processor implements the method described above when executing the computer program.

[0045] Based on the same inventive concept, this disclosure also provides a non-transitory computer-readable storage medium that stores computer instructions for causing a computer to perform the method described above.

[0046] As can be seen from the above, the method provided in this application combines high-precision manual monitoring and satellite positioning-acquired second coordinate data with first coordinate data based on a digital surface model to correct the raster data, effectively compensating for the insufficient accuracy of a single data source. Compared to traditional methods that rely solely on UAV monitoring, this method can more accurately reflect subtle deformations of the mining area surface, greatly improving the accuracy of monitoring data and providing reliable data support for safe production and disaster early warning in mining areas. Compared to relying solely on high-precision manual monitoring and satellite positioning monitoring, this method only requires acquiring a few high-precision points, reducing the scope and frequency of high-precision measurements, thereby lowering the overall monitoring cost.

[0047] Furthermore, addressing the issue of significant differences in dimensionality and spatiotemporal resolution in traditional multi-source data, the rasterization of digital surface models facilitates the unification of discrete second coordinate data (high-precision coordinate data) into a raster coordinate system, eliminating fusion barriers caused by geographic spatial heterogeneity. Additionally, since the raster data and the second coordinate data of monitoring points can originate from the same point in time, it effectively avoids the "time misalignment" problem caused by different collection times in traditional multi-source data. Moreover, rasterization facilitates the correction of raster data through the model, enabling both fine-tuning of local errors and rapid updates of the entire data area through algorithmic iteration. Attached Figure Description

[0048] To more clearly illustrate the technical solutions in this application or related technologies, the drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0049] Figure 1 A flowchart of the data fusion method is shown for an embodiment of this application;

[0050] Figure 2 This application provides a map showing the location of the monitoring points in an embodiment.

[0051] Figure 3 This is a schematic diagram showing the grid corresponding to monitoring point A2 in an embodiment of this application;

[0052] Figure 4 This application provides an embodiment showing the elevation map corresponding to monitoring point A2;

[0053] Figure 5 This is a schematic diagram of the hardware structure of an electronic device according to an embodiment of this application. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with specific embodiments and the accompanying drawings.

[0055] It should be noted that, unless otherwise defined, the technical or scientific terms used in the embodiments of this application should have the ordinary meaning understood by one of ordinary skill in the art to which this application pertains. The terms "first," "second," and similar terms used in the embodiments of this application do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are only used to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0056] As mentioned in the background, surface subsidence monitoring plays a crucial role in the field of geological safety, as its results directly impact infrastructure safety, ecological environment protection, and the safety of people's lives and property. Traditional monitoring methods, such as leveling, InSAR (Inductive Aperture Radar Interferometry), and GPS (Global Positioning System), each have significant limitations. While leveling offers high accuracy, it requires substantial manpower and has a long measurement cycle, making real-time dynamic monitoring difficult. InSAR technology, although capable of acquiring surface deformation information over a wide area, is affected by atmospheric effects and vegetation cover, significantly reducing its accuracy in complex terrain areas. GPS provides precise positioning but can only acquire displacement data at discrete points, failing to fully reflect the continuous deformation characteristics of the surface. The shortcomings of these individual monitoring methods in terms of data accuracy and application scope make it difficult to comprehensively capture the dynamic changes in surface subsidence and meet the ever-increasing demand for high-precision, real-time monitoring.

[0057] Given the limitations of traditional monitoring methods, multi-source data fusion technology has emerged and is gradually becoming an important development direction for land subsidence monitoring. This technology integrates data acquired through different monitoring methods, aiming to leverage the strengths of each method and compensate for their respective weaknesses, thereby achieving more comprehensive and accurate monitoring of land subsidence. However, in practical applications, multi-source data fusion technology faces many severe challenges.

[0058] The primary challenge lies in the significant differences and geospatial heterogeneity of data acquired through different monitoring methods. For example, optical remote sensing imagery provides high-resolution visual information about the Earth's surface, but its temporal resolution is relatively low; while InSAR data can monitor periodic deformation, its spatial resolution is relatively limited. Furthermore, the measurement principles, coordinate systems, and time bases of different sensors vary, resulting in vast differences in dimensionality and spatiotemporal resolution, making direct assimilation and unification difficult. This data heterogeneity not only increases the difficulty of data preprocessing but may also lead to information loss or erroneous fusion, severely impacting the accuracy and reliability of monitoring results.

[0059] Another key issue is how to choose a suitable method to achieve efficient and accurate multi-source information fusion. Since the data acquired by each monitoring method has different characteristics, the applicable fusion methods also differ. For example, for complementary point data (such as GPS) and regional data (such as InSAR), how to organically combine the two, retaining the high-precision advantage of point data while utilizing the spatial coverage characteristics of regional data, is a complex technical challenge.

[0060] The following is in conjunction with the appendix Figures 1-5 The embodiments of this application will be described in detail below.

[0061] like Figure 1 As shown, a method for fusing multi-source monitoring data of surface subsidence in mining areas includes the following steps:

[0062] S100: Acquire monitoring image data of the preset research area by the UAV, and generate a digital surface model based on the monitoring image data;

[0063] Specifically, within the pre-defined study area, drone flight routes are planned based on the mine's topography, mining scope, and monitoring needs. Using drones equipped with high-resolution optical cameras or lidar, multiple sorties are conducted along the set routes at multiple angles to acquire high-definition monitoring image data covering the entire study area. The acquired image data is then processed using photogrammetry and other processing techniques to generate a high-precision digital surface model (DSM), which visually reflects the three-dimensional morphological characteristics of the study area's surface.

[0064] S200: Cut the digital surface model into multiple grids, and determine the grid data of each grid based on the digital surface model;

[0065] Specifically, based on the generated digital surface model, and according to the monitoring accuracy requirements and data processing capabilities, the model is divided into multiple grids of appropriate size. Each grid acts as a data unit; by performing gridding processing on the digital surface model, the grid data contained within each grid is determined. This data records the spatial location information of the land surface within the grid area, providing basic data for subsequent monitoring and analysis.

[0066] S300: Detection points are deployed in at least a portion of the grid, and first coordinate data and second coordinate data of the detection points are acquired;

[0067] Specifically, monitoring points are scientifically and rationally deployed in at least some key grids, taking into account factors such as the impact range of mining activities and geological structural characteristics. For these monitoring points, on the one hand, their first coordinate data is determined based on a digital surface model. This first coordinate data is a specific data point within the corresponding grid, providing approximate location information of the monitoring point within the model. On the other hand, high-precision manual monitoring methods are employed, such as leveling and angle measurements using a total station, or second coordinate data of the monitoring points obtained using a high-precision satellite positioning system (such as RTK-GPS). Because the second coordinate data utilizes more precise measurement methods and equipment, its accuracy is significantly higher than the first coordinate data, more accurately reflecting the actual spatial location of the monitoring point.

[0068] S400: Based on the first coordinate data and the second coordinate data of the monitoring point, the grid data of each grid is corrected to obtain the corrected grid data;

[0069] Specifically, by comparing and analyzing the first and second coordinate data of the same monitoring point, the coordinate error value is calculated. These error values ​​are then applied to the raster data of each grid cell, and each grid cell is corrected one by one. For example, for each grid cell, its coordinates are adjusted according to its spatial relationship with the monitoring point and the error model, thereby obtaining the corrected grid data, making the corrected data closer to the actual surface morphology.

[0070] S500: Identify all corrected raster data as multi-source monitoring and fusion data of surface subsidence in the mining area;

[0071] The first coordinate data is determined based on the digital surface model, while the second coordinate data is not determined based on the digital surface model, and the accuracy of the second coordinate data is higher than that of the first coordinate data.

[0072] In addition, the second coordinate data includes coordinate data obtained from manual monitoring and coordinate data obtained from satellite positioning.

[0073] In this embodiment, by combining the second coordinate data obtained through high-precision manual monitoring and satellite positioning with the first coordinate data based on a digital surface model, the raster data is corrected, effectively compensating for the insufficient accuracy of a single data source. Compared to traditional methods that rely solely on UAV monitoring, this method can more accurately reflect subtle deformations of the mining area's surface, greatly improving the accuracy of monitoring data and providing reliable data support for safe production and disaster early warning in the mining area. Compared to relying solely on high-precision manual monitoring and satellite positioning monitoring, this method only requires acquiring a few high-precision points, reducing the scope and frequency of high-precision measurements, thereby lowering the overall monitoring cost.

[0074] Furthermore, addressing the issue of significant differences in dimensionality and spatiotemporal resolution in traditional multi-source data, the rasterization of digital surface models facilitates the unification of discrete second coordinate data (high-precision coordinate data) into a raster coordinate system, eliminating fusion barriers caused by geographic spatial heterogeneity. Additionally, since the raster data and the second coordinate data of monitoring points can originate from the same point in time, it effectively avoids the "time misalignment" problem caused by different collection times in traditional multi-source data. Moreover, rasterization facilitates the correction of raster data through the model, enabling both fine-tuning of local errors and rapid updates of the entire data area through algorithmic iteration.

[0075] In some embodiments, in step S400, the process of correcting the raster data of each grid cell based on the first coordinate data and the second coordinate data of the monitoring point to obtain corrected raster data includes:

[0076] S410: In the pre-constructed improved radial basis function interpolation algorithm model, the raster data of each grid is corrected based on the first coordinate data and the second coordinate data of the monitoring point to obtain the corrected raster data.

[0077] In some embodiments, the pre-construction process of the improved radial basis function interpolation algorithm model in step S410 includes the following steps:

[0078] S411: Determine the terrain roughness and tilt of each grid cell;

[0079] Specifically, based on the raster data of the Digital Surface Model (DSM), the terrain roughness and tilt of each raster are calculated using the local neighborhood analysis method.

[0080] S412: Determine the mean terrain roughness of the study area based on the terrain roughness of all the grids;

[0081] Specifically, the terrain roughness values ​​of all grid cells are aggregated, and the arithmetic mean of the terrain roughness of the entire study area is calculated. This mean serves as a benchmark indicator of the regional terrain characteristics and is used for subsequent algorithm weight allocation.

[0082] S413: Perform the following steps for each of the grid cells:

[0083] S413a: Determine the shape parameters of a grid based on its terrain roughness and tilt.

[0084] Specifically, the shape parameters of the grid are determined by combining the terrain roughness and tilt of the grid, and these shape parameters affect the range of the subsequent Gaussian function.

[0085] S413b: Based on the shape parameters of the grid, the grid data of the grid and the first coordinate data of the corresponding target monitoring point, a Gaussian function model is constructed; each grid corresponds to at least one target monitoring point, the target monitoring point being a monitoring point located within the grid or a monitoring point located around the grid;

[0086] S413c: Based on the grid data of this grid and the first coordinate data of the corresponding target monitoring point, construct a multiple quadratic function model;

[0087] Specifically, the Gaussian function is mainly used to describe the degree of influence of monitoring points on each point within the grid. It is suitable for flat areas. In mining area monitoring scenarios, it can reflect the impact of monitoring points on the surrounding areas from a macro perspective. For example, in areas with relatively simple terrain and sparse distribution of monitoring points, the Gaussian function can effectively utilize limited monitoring point information to perform preliminary correction and trend prediction of grid data.

[0088] The quadratic function model is mainly used to describe the influence of monitoring points on each point within the grid. It is suitable for gentle slopes or gully areas. In areas with complex terrain and obvious deformation characteristics, such as the area around mining subsidence areas, the quadratic function can capture subtle changes and local anomalies in the terrain more precisely, and make more accurate corrections to the grid data, thus making up for the limitations of the Gaussian function in characterizing the details of complex terrain.

[0089] S413d: Based on the mean terrain roughness and the terrain roughness of the grid, determine the first weight for the Gaussian function model and the second weight for the multiple quadratic function model;

[0090] S413e: Determine the influence factor of the monitoring point on the grid based on the first weight, Gaussian function model, second weight and multiple quadratic function model;

[0091] Specifically, the sum of the first and second weights is 1. Weights are assigned based on the mean terrain roughness of the study area and the terrain roughness of the current raster. For raster with flat terrain, the first weight of the Gaussian function is higher; for raster with complex terrain, the second weight of the quadratic function is increased to enhance the ability to correct local details. An influence factor is obtained through the first weight, Gaussian function model, second weight, and quadratic function model. This influence factor is used to quantify the degree of influence of different monitoring points on raster correction.

[0092] S413f: Determine the elevation residual of the first and second coordinate data of the monitoring point corresponding to at least one target monitoring point;

[0093] S413g: Based on the first coordinate data of at least one target monitoring point, the elevation residual, and the influence factor on the grid, an improved radial basis function interpolation algorithm model is constructed.

[0094] Specifically, in this improved radial basis function interpolation algorithm model, the raster data of each grid is corrected based on the first and second coordinate data of the monitoring point to obtain the corrected raster data.

[0095] In this embodiment, the traditional radial basis function (RBF) interpolation algorithm estimates unknown point data by constructing a linear combination of basis functions. However, when processing multi-source heterogeneous monitoring data in mining areas, its interpolation accuracy often fails to meet practical requirements due to the complex and varied terrain and diverse data sources. The improved RBF interpolation algorithm model optimizes the data correction process from the root by introducing an adaptive weight allocation mechanism and terrain constraints.

[0096] This model dynamically adjusts the weighting ratio of Gaussian and quadratic functions based on terrain features to accurately calculate the influence factor of monitoring points on the raster. Specifically, the simpler the terrain, the higher the weight of the Gaussian function (emphasizing the global picture); the more complex the terrain, the higher the weight of the quadratic function (strengthening the local picture), achieving a dynamic balance between "global trends and local details." In simple terrain areas, the influence factor integrates more global information from monitoring points, making data correction more consistent with macro-terrain trends; in complex terrain areas, the influence factor focuses on local details, improving the correction accuracy for minor terrain changes. When actually correcting raster data, the algorithm comprehensively considers multiple factors such as the spatial distance between raster points and monitoring points, terrain features, and influence factors, dynamically adjusting the interpolation weights to ensure that the coordinates of each raster point are optimally corrected. This improved RBF interpolation algorithm model can stably control the correction error of raster data to the millimeter level, effectively solving the accuracy bottleneck of traditional algorithms in complex terrain, providing high-precision and highly adaptable data support for monitoring surface subsidence in mining areas, and greatly improving the reliability and practical value of monitoring data.

[0097] In some embodiments, in step S411, determining the terrain roughness and tilt of each of the grid cells includes:

[0098] S411a: Obtain the raster data of a grid cell and all its adjacent grid cells, and based on the raster data of the grid cell and all its adjacent grid cells, obtain the elevation data of the grid cell and all its adjacent grid cells.

[0099] Specifically, the elevation data of a raster can be obtained by summing the elevation values ​​of all data within the raster and then dividing by the total number of data points. This method is simple and direct, and can effectively reflect the overall elevation level within the raster area. Alternatively, mathematical functions can be used to fit the data within the raster, such as polynomial fitting or surface fitting, to construct a mathematical model that describes the terrain of the raster. The elevation values ​​representing the raster are then extracted from this model. This method is suitable for areas with complex terrain and irregular data distribution, and can more accurately reflect the terrain change trends within the raster, providing more accurate basic data for subsequent terrain analysis and data processing.

[0100] S411b: Obtain the average elevation data based on the elevation data of this grid and all adjacent grids;

[0101] Specifically, all adjacent rasters of the target raster are defined as all rasters within a 3×3 window area, extending one row and one column horizontally and vertically from the target raster itself as the center. This window encompasses the target raster itself and its eight directly adjacent raster units (including orthogonal and diagonal adjacent rasters), collectively forming the basic data unit for local terrain feature analysis.

[0102] S411c: Determine the terrain roughness of the grid based on the elevation data and average elevation data of the grid and all adjacent grids.

[0103] Specifically, the terrain roughness of the raster is expressed by the formula:

[0104] ;

[0105] in, For the first The terrain roughness of the grid; For the first Within a 3x3 window centered on the grid, the row and column offsets are: adjacent grid elevations ; The average elevation value within a 3×3 window.

[0106] In some embodiments, in step S411, determining the terrain roughness and tilt of each of the grid cells further includes:

[0107] S411d: Obtain the raster data of a grid cell and its neighboring grid cells in the x-direction; based on the raster data of the grid cell and its neighboring grid cells in the x-direction, obtain the elevation data of the grid cell and its neighboring grid cells in the x-direction.

[0108] S411e: Based on the elevation data of the grid and its adjacent grids in the x-direction, obtain the height gradient of the grid in the x-direction;

[0109] S411f: Obtain the raster data of a grid cell and its neighboring grid cells in the y-direction; based on the raster data of the grid cell and its neighboring grid cells in the y-direction, obtain the elevation data of the grid cell and its neighboring grid cells in the y-direction.

[0110] S411g: Based on the elevation data of the grid and its adjacent grids in the y direction, obtain the height gradient of the grid in the y direction;

[0111] S411h: Obtain the elevation scaling factor of the digital surface model;

[0112] S411i: Determine the tilt of the grid based on the height gradient of the grid in the x-direction and the height gradient in the y-direction, as well as the elevation scaling factor.

[0113] The x-direction and y-direction are located on the same horizontal plane.

[0114] Specifically, the grid tilt is expressed by the formula:

[0115] ;

[0116] ;

[0117] ;

[0118] in, For the first The tilt of the grid; This represents the height gradient of the grid in the x-direction; For this grid, the height gradient in the y-direction; This is the elevation scaling factor (DSM elevation units are in meters, therefore...). =1); and f x Related h i+1，j and h i-1，j For the first Elevation data of two adjacent rasters in the x-direction; and f y Related h i , j+1 and hi，j-1 For the first Elevation data of two adjacent grid cells in the y-direction; resolution is the same as that of the UAV monitoring image.

[0119] In some embodiments, determining the mean terrain roughness of the study area based on the terrain roughness of all the grids includes:

[0120] Based on the terrain roughness of all the grids and the mean terrain roughness of the study area, the standard deviation of the terrain roughness of the study area is determined;

[0121] The method of determining the shape parameters of a grid based on its terrain roughness and tilt also includes:

[0122] The shape parameters of the grid are determined based on the mean and standard deviation of the terrain roughness of the study area, as well as the terrain roughness and tilt of the grid.

[0123] Specifically, the shape parameters of the grid are expressed by the formula:

[0124] ;

[0125] in, For the first The shape parameters of the grid; The mean topographic roughness of the study area; This represents the standard deviation of terrain roughness.

[0126] In this embodiment, the shape parameter directly determines the influence range and decay rate of the Gaussian function. This embodiment determines the shape parameter of the Gaussian function by studying the standard deviation and mean of terrain roughness in the research area, making the influence range of the Gaussian function more closely match the actual terrain. In areas with complex terrain and large standard deviation of roughness, the shape parameter obtained by the relevant calculation is smaller, resulting in faster decay of the Gaussian function and a more concentrated influence range. This allows for precise focus on local complex terrain changes and capture of subtle terrain features. When calculating the influence factor of monitoring points on the target grid, this concentrated influence range will prioritize integrating data information from nearby monitoring points around the target grid, ensuring that the obtained influence factor more accurately reflects the local correlation and change patterns of monitoring data under complex terrain. In areas with simple terrain and small standard deviation of roughness, the shape parameter is larger, the Gaussian function decays more slowly, and the influence range is wider, focusing on the control of macro terrain trends. At this time, the influence range of the Gaussian function is large, and when calculating the influence factor, more emphasis is placed on using data from monitoring points within a larger range around the target grid. The Gaussian function model enhances the ability to characterize the overall terrain trend, while combining multiple quadratic function models to assist in optimizing local details.

[0127] In some embodiments, the determination of the influence factor of the monitoring point on the grid based on the first weight, Gaussian function model, second weight, and multiple quadratic function model is expressed by the following formula:

[0128] ;

[0129] in, The influence factor of the i-th target monitoring point on the grid; It is a Gaussian function model. These are the shape parameters of the grille; It is a model of multiple quadratic functions; The Euclidean distance between the grid and the corresponding target monitoring point on the horizontal plane is calculated using the grid data of the grid and the first coordinate data of the corresponding target monitoring point. It is the first weight; is the second weight, and c is the rate constant.

[0130] Specifically, Gaussian function model It is constructed based on the shape parameters of the grid, the grid data of the grid, and the first coordinate data of the corresponding target monitoring point, wherein, This refers to the Euclidean distance on the horizontal plane between the grid and the corresponding target monitoring point, calculated using the grid data and the first coordinate data of the target monitoring point. Further, the x and y coordinates of all data within the grid are averaged to obtain the coordinates of the grid center point. The Euclidean distance on the horizontal plane between this grid and the corresponding target monitoring point is the distance between the grid center point and the target monitoring point on the horizontal plane. i Expressed using a formula: Where (x, y) are the coordinates of the center point of the grid, (x...y...) i y i ) represents the coordinates of the i-th target monitoring point. The model is a quadratic function model, where c is a fixed parameter of the quadratic function that controls the sensitivity to nearby abrupt changes; for example, c=5. The first weight for the Gaussian function model and the second weight for the quadratic function model are based on the mean terrain roughness and the terrain roughness of the raster, and are expressed by the following formula: .

[0131] In this embodiment, weights are dynamically assigned based on the mean terrain roughness and the roughness of the raster itself. The simpler the terrain, the higher the weight of the Gaussian function (emphasizing global factors); the more complex the terrain, the higher the weight of the quadratic function (strengthening local factors), achieving a dynamic balance between "global trends and local details." Through dynamic weight allocation, the influence factor can more accurately reflect the impact of terrain features on the monitoring data. In simple terrain areas, the influence factor integrates more global information from monitoring points, making data correction more consistent with macro-terrain trends; in complex terrain areas, the influence factor focuses on local details, improving the accuracy of corrections for minor terrain changes.

[0132] In some embodiments, an improved radial basis function interpolation algorithm model is constructed based on the first coordinate data of at least one target monitoring point, the elevation residual, and the influence factor on the grid.

[0133] Specifically, based on the first coordinate data of at least one target monitoring point, the elevation residual, and the influence factors on the raster, an improved radial basis function interpolation algorithm model is constructed, including:

[0134] Based on the first coordinate data, elevation residual, and influence factor of at least one target monitoring point, determine the contribution weight and polynomial coefficient of each monitoring point to the grid.

[0135] Based on the influence factors, contribution weights, and polynomial coefficients of each monitoring point on the raster, the corrected elevation residual of the raster is obtained.

[0136] Based on the elevation value of the grid and the corrected elevation residual, an improved radial basis function interpolation algorithm model is constructed.

[0137] The contribution weights and polynomial coefficients of this grid are solved using the following set of constraint equations:

[0138] ;

[0139] in, The influence factor of the i-th monitoring point on the target grid is a combination of Gaussian and quadratic function characteristics.

[0140] di: the horizontal Euclidean distance between the i-th monitoring point and the target grid; Δh i The elevation residual of the i-th grid cell (the deviation between the actual elevation and the initial model); (x j y j ): The planar coordinates of the j-th monitoring point; N: The total number of monitoring points; a0, a1, a2: Polynomial coefficients (a = [a0, a1, a2]) T ); λ j : The contribution weight of the j-th monitoring point (λ=[λ1, λ2, ..., λj) N ] T ); Let be the planar coordinates of any grid.

[0141] To facilitate numerical solutions, the constraint equations are transformed into block matrix form:

[0142] ;

[0143] Where λ is the weight vector; a is the polynomial coefficient vector; Φ is an M×N matrix (M is the total number of grids in the study area, N is the number of monitoring points); P is an N×3 matrix, each row corresponding to a monitoring point, with content [1, x... j y j ]; It is an M×1 vector that stores the elevation residuals of each grid cell (Δh=[Δh1,Δh2,…,Δh...). M ] T );

[0144] Based on the elevation value of the raster and the corrected elevation residual, an improved radial basis function interpolation algorithm model is constructed, which is expressed as follows:

[0145] ;

[0146] ;

[0147] For each grid cell, the correction amount is calculated using k high-precision points, and the correction is performed point by point to calculate the corrected elevation residual. Then, the corrected elevation residuals and the raster elevations are compared. Add them together to get the fusion value. Generate a fused high-precision elevation, with the grid coordinates as follows: , This is the merged elevation value.

[0148] In this embodiment, the improved radial basis function interpolation algorithm model addresses the underfitting problem of traditional algorithm models in complex terrains of mining areas (such as goafs, steep slopes, and gullies). By fusing the wide-area trend capture of Gaussian functions with the local detail correction of multiple quadratic functions, the terrain of the mining area is accurately reconstructed.

[0149] The above embodiments will be described below with reference to specific examples:

[0150] Example 1

[0151] To verify the feasibility of the method, the 12308 working face of Haohua Hongqingliang Mining Co., Ltd. in Ordos City, Inner Mongolia, was selected as the research object. Multi-source monitoring data collection was carried out on June 1, 2023, and data from UAV aerial survey, GNSS positioning and manual monitoring were obtained.

[0152] The DSM digital surface model of the mining area was obtained using UAV photogrammetry technology. The research area was a 250×250m area.

[0153] The UAV DSM digital surface model was divided into 158,840 square grids of 0.65 × 0.65 m each. Data from all grids was extracted. Ten monitoring points were set up within each grid. The second coordinate data for six monitoring points (A1-A6) came from GNSS monitoring, and the second coordinate data for four monitoring points (B1-B4) came from manual monitoring. All data were uniformly processed, with the plane coordinate system consistent with the 1954 Beijing Coordinate System, and the elevation datum uniformly adopted from the 1985 National Elevation Datum. The UAV imagery, GNSS data, and the locations of the manual monitoring points are shown below. Figure 2 As shown in Table 1, the second coordinate data corresponding to GNSS and manual monitoring points are as follows.

[0154] Table 1 Coordinates of GNSS and Manual Measurement Points

[0155]

[0156] Taking GNSS monitoring point A2 as an example, such as Figure 3 and Figure 4 As shown in Table 2, the UAV grid number corresponding to GNSS monitoring point A2 is 54916, the elevation of A2 is 1424.34m, the terrain roughness index TRI is 0.12, and the tilt is... .

[0157] To verify the accuracy of the fused data obtained by this method and its applicability under different terrains, monitoring points A2 in grassland terrain and A4 in gully terrain were excluded from fusion. Based on the first and second coordinate data of eight monitoring points (GNSS monitoring points A1, A3, A5, and A6, and manual monitoring points B1, B2, B3, and B4), a modified radial basis function interpolation algorithm was used to correct the global raster data, resulting in fused data, including the fused data (corrected first coordinate data) for monitoring points A2 and A4. Finally, the second coordinate data of monitoring points A2 and A4 were compared with the fused data (corrected first coordinate data) to verify the accuracy of the fused data under different terrains. The second coordinate data, first coordinate data, corrected first coordinate data, and absolute value of accuracy improvement for grassland (A2) and gully terrain (A4) are shown in the table below.

[0158] Table 2 Statistical Table of Fusion Results

[0159]

[0160] As shown in Table 2, due to the noise of the drone itself, the elevation difference between the first and second coordinate data of monitoring point A2 in the grassland area was relatively large, at 3.24m. However, after correction using this method, the difference between the elevation values ​​of the first and second coordinate data was 0.984m, which is an improvement of about 69.6% compared to the elevation value of the first coordinate data. In the gully area, the elevation difference between the first and second coordinate data of monitoring point A4 was 0.220m due to the influence of the gully terrain. After correction using this method, the difference between the elevation values ​​of the first and second coordinate data was 0.069m, which is an improvement of about 68%.

[0161] In summary, the fusion data obtained by the multi-source monitoring data fusion method for surface subsidence in mining areas based on the improved RBF proposed in this invention has significantly improved the accuracy of the fused data compared with the original data from UAVs. This proves the feasibility of the multi-source data fusion method proposed in this invention in terms of accuracy, and it is still applicable to grassland and gully landforms.

[0162] It should be noted that the method in this embodiment can be executed by a single device, such as a computer or server. The method can also be applied in a distributed scenario, where multiple devices cooperate to complete the task. In such a distributed scenario, one of these devices may execute only one or more steps of the method in this embodiment, and the multiple devices will interact with each other to complete the method described.

[0163] It should be noted that the above description describes some embodiments of this application. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in a different order than that shown in the above embodiments and still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0164] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this application also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the data fusion method described in any of the above embodiments.

[0165] Figure 5This embodiment illustrates a more specific hardware structure of an electronic device. The device may include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040, and a bus 1050. The processor 1010, memory 1020, input / output interface 1030, and communication interface 1040 are interconnected internally via the bus 1050.

[0166] The processor 1010 can be implemented using a general-purpose CPU (Central Processing Unit), microprocessor, application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.

[0167] The memory 1020 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 can store the operating system and other applications. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program code is stored in the memory 1020 and is called and executed by the processor 1010.

[0168] The input / output interface 1030 is used to connect input / output modules to realize information input and output. Input / output modules can be configured as components within the device (not shown in the figure) or externally connected to the device to provide corresponding functions. Input devices may include keyboards, mice, touchscreens, microphones, various sensors, etc., while output devices may include displays, speakers, vibrators, indicator lights, etc.

[0169] The communication interface 1040 is used to connect a communication module (not shown in the figure) to enable communication between this device and other devices. The communication module can communicate via wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.).

[0170] Bus 1050 includes a pathway for transmitting information between various components of the device, such as processor 1010, memory 1020, input / output interface 1030, and communication interface 1040.

[0171] It should be noted that although the above-described device only shows the processor 1010, memory 1020, input / output interface 1030, communication interface 1040, and bus 1050, in specific implementations, the device may also include other components necessary for normal operation. Furthermore, those skilled in the art will understand that the above-described device may only include the components necessary for implementing the embodiments of this specification, and not necessarily all the components shown in the figures.

[0172] The electronic devices described above are used to implement the corresponding data fusion methods in any of the foregoing embodiments and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0173] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this application also provides a non-transitory computer-readable storage medium storing computer instructions for causing the computer to execute the data fusion method as described in any of the above embodiments.

[0174] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.

[0175] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of this application is limited to these examples; under the concept of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of the embodiments of this application as described above, which are not provided in detail for the sake of brevity.

[0176] Additionally, to simplify the description and discussion, and to avoid obscuring the embodiments of this application, the well-known power / ground connections to integrated circuit (IC) chips and other components may or may not be shown in the provided drawings. Furthermore, the apparatus may be shown in block diagram form to avoid obscuring the embodiments of this application, and this also takes into account the fact that the details of the implementation of these block diagram apparatuses are highly dependent on the platform on which the embodiments of this application will be implemented (i.e., these details should be fully understood by those skilled in the art). While specific details (e.g., circuits) have been set forth to describe exemplary embodiments of this application, it will be apparent to those skilled in the art that the embodiments of this application can be implemented without these specific details or with variations thereof. Therefore, these descriptions should be considered illustrative rather than restrictive.

[0177] Although this application has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed.

[0178] The embodiments of this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of the claims of this application. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the embodiments of this application should be included within the protection scope of this application.

Claims

1. A method for fusing multi-source monitoring data of surface subsidence in mining areas, characterized in that, include: Acquire monitoring image data of a pre-defined research area by a drone, and generate a digital surface model based on the monitoring image data; The digital surface model is cut into multiple grids, and the grid data of each grid is determined based on the digital surface model; Monitoring points are deployed in at least a portion of the grid, and the first and second coordinate data of the monitoring points are acquired. In the pre-constructed improved radial basis function interpolation algorithm model, the raster data of each grid is corrected based on the first and second coordinate data of the monitoring point to obtain the corrected raster data. The corrected raster data of all graticules were identified as multi-source monitoring and fusion data of surface subsidence in the mining area; The first coordinate data is determined based on the digital surface model, while the second coordinate data is not determined based on the digital surface model. The accuracy of the second coordinate data is higher than that of the first coordinate data. The pre-construction process of the improved radial basis function interpolation algorithm model includes: Determine the terrain roughness and tilt of each grid cell; The mean topographic roughness of the study area is determined based on the topographic roughness of all grids. Perform the following steps for each grid cell: Determine the shape parameters of a grid based on its terrain roughness and tilt. Based on the shape parameters of the grid, the grid data of the grid, and the first coordinate data of the corresponding target monitoring point, a Gaussian function model is constructed; each grid corresponds to at least one target monitoring point, which is either a monitoring point located within the grid or a monitoring point located around the grid; Based on the grid data of the grid and the first coordinate data of the corresponding target monitoring point, a multi-quadratic function model is constructed; Based on the mean terrain roughness and the terrain roughness of the grid, determine the first weight for the Gaussian function model and the second weight for the multiple quadratic function model; Based on the first weight, Gaussian function model, second weight and multiple quadratic function model, the influence factor of the monitoring point on the grid is determined; Determine the elevation residuals of the first and second coordinate data of at least one target monitoring point; An improved radial basis function interpolation algorithm model is constructed based on the first coordinate data of at least one target monitoring point, the elevation residual, and the influence factor on the grid.

2. The method according to claim 1, characterized in that, The second coordinate data includes coordinate data obtained through manual monitoring and coordinate data obtained through satellite positioning.

3. The method according to claim 1, characterized in that, The determination of the terrain roughness and tilt of each grid cell includes: Obtain the raster data of a given grid cell and all its adjacent grid cells. Based on the raster data of the given grid cell and all its adjacent grid cells, obtain the elevation data of the given grid cell and all its adjacent grid cells. Based on the elevation data of this grid and all adjacent grids, obtain the average elevation data; The terrain roughness of the grid is determined based on the elevation data and average elevation data of the grid and all adjacent grids.

4. The method according to claim 1, characterized in that, The determination of the terrain roughness and tilt of each grid cell also includes: Obtain raster data of a grid cell and its neighboring grid cells in the x-direction. Based on the raster data of the grid cell and its neighboring grid cells in the x-direction, obtain the elevation data of the grid cell and its neighboring grid cells in the x-direction. Based on the elevation data of the grid and its adjacent grids in the x-direction, the height gradient of the grid in the x-direction is obtained; Obtain raster data of a grid cell and its neighboring grid cells in the y-direction. Based on the raster data of the grid cell and its neighboring grid cells in the y-direction, obtain the elevation data of the grid cell and its neighboring grid cells in the y-direction. Based on the elevation data of the grid and its adjacent grids in the y-direction, the height gradient of the grid in the y-direction is obtained; Obtain the elevation scaling factor of the digital surface model; The tilt of the grid is determined based on the height gradient of the grid in the x-direction and the height gradient in the y-direction, as well as the elevation scaling factor. The x-direction and y-direction are located on the same horizontal plane.

5. The method according to claim 1, characterized in that, The determination of the mean terrain roughness of the study area based on the terrain roughness of all grids includes: Based on the terrain roughness of all grids and the mean terrain roughness of the study area, the standard deviation of terrain roughness of the study area is determined. The method of determining the shape parameters of a grid based on its terrain roughness and tilt also includes: The shape parameters of the grid are determined based on the mean and standard deviation of the terrain roughness of the study area, as well as the terrain roughness and tilt of the grid.

6. The method according to claim 1, characterized in that, The influence factor of the monitoring point on the grid is determined based on the first weight, Gaussian function model, second weight, and multiple quadratic function model, and is expressed by the following formula: Wherein, φMixed(d i ) represents the influence factor of the i-th target monitoring point on the grid; It is a Gaussian function model; It is a model of multiple quadratic functions; d i The distance between the grid and the corresponding target monitoring point on the horizontal plane is calculated using the grid data of the grid and the first coordinate data of the corresponding target monitoring point; α is the first weight; 1-α is the second weight; and c is the rate constant.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the program, it implements the method as described in any one of claims 1 to 6.

8. A non-transitory computer-readable storage medium storing computer instructions, characterized in that, The computer instructions are used to cause the computer to perform the method according to any one of claims 1 to 6.

Images