Mining area settlement monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion

Abstract

The invention provides a mining area settlement monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion, and belongs to the technical field of mining area settlement monitoring and grouting treatment evaluation. In order to solve the technical problems of incomplete monitoring dimension, disjunction of monitoring and treatment, and insufficient timeliness and automation degree in the prior art, monitoring calculation equipment is used for obtaining time sequence SAR image data covering a target mining area, thermal noise removal, radiometric calibration and multi-view processing are performed on the SAR image data, and then an SBAS-InSAR technology is adopted for processing, so that the time sequence SAR image data is obtained. Obtaining a deformation rate field and an accumulated deformation quantity of the whole mining area, identifying a significant settlement area and a decoherence area, and generating a superposed vector diagram of the key area and the decoherence area; the method comprises the following steps of: carrying laser radar equipment by an unmanned aerial vehicle, and respectively executing an aerial survey task of the unmanned aerial vehicle at two time points of an initial settlement period and a settlement development period to obtain vertical settlement volumes of a key area and an incoherent area; the method is applied to mining area settlement monitoring and treatment evaluation.

Description

TECHNICAL FIELD

[0001] The application provides a mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion, and belongs to the technical field of mine area subsidence monitoring and grouting treatment evaluation. BACKGROUND

[0002] The subsidence of a mined-out area in a mine area is a major hidden danger threatening the safety production of the mine and the stability of the surrounding geological environment. At present, off-layer grouting is used as an effective treatment method for the subsidence of the mined-out area, and the treatment effect depends to a great extent on the accurate grasp of the subsidence range, deformation characteristics and development trend. However, there are still the following technical bottlenecks in the monitoring of the subsidence of the mined-out area in the mine area and the grouting treatment.

[0003] Incomplete monitoring dimension: traditional single monitoring technology cannot obtain complete deformation field information. Although discrete point measurement (such as GNSS and leveling) has high precision, it cannot realize full-field coverage. Satellite remote sensing technology can monitor a large range, but it is easy to appear incoherence due to the influence of vegetation coverage and mining activities in the mine area, and can only provide one-dimensional deformation along the radar line-of-sight direction, and cannot accurately reflect the true three-dimensional deformation characteristics of the ground.

[0004] Monitoring and treatment disconnection: the existing technical system divides deformation monitoring and grouting treatment into two relatively independent links. The monitoring data often stays at the deformation description level, and cannot be directly and quantitatively associated with the grouting engineering design parameters. The determination of the grouting target area depends on engineering experience, lacks scientific data support, and leads to a certain blindness in the treatment engineering, resulting in poor treatment effect.

[0005] Insufficient timeliness and automation: the traditional monitoring method has long cycle and low efficiency, and it is difficult to timely reflect the rapid deformation characteristics of the mine area. The data processing and analysis process relies on manual intervention, and the automation degree is not high, which cannot meet the needs of real-time early warning and rapid response of the safety production of the mine area. SUMMARY

[0006] In order to solve the technical problems in the background art, the technical scheme adopted by the application is as follows: a mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion is provided, which comprises the following monitoring and evaluation steps:

[0007] Step S1: using a monitoring and calculation device to obtain time-series SAR image data covering a target mine area, preferentially selecting C-band Sentinel-1 data or X-band TerraSAR-X data, performing thermal noise removal, radiation calibration and multi-view processing on the SAR image data, and then using SBAS-InSAR technology for processing to obtain InSAR deformation field data of the entire mine area, including deformation rate field and cumulative deformation, identifying significant subsidence areas and incoherent areas, and generating superimposed vector maps of key areas and incoherent areas;

[0008] Step S2: Using the superimposed vector diagram of the key area and the incoherent area generated in step S1, the monitoring computing device automatically plans the aerial survey range and route of the unmanned aerial vehicle LiDAR, and the unmanned aerial vehicle carries the laser radar device to perform the unmanned aerial vehicle LiDAR aerial survey task at two time points of the initial settlement T1 and the settlement development T2, respectively, to obtain the vertical settlement Dz of the key area and the incoherent area;

[0009] Step S3: Using the monitoring computing device, the InSAR deformation field data obtained in step S1 is fused with the vertical settlement Dz obtained in step S2 to calculate the east-west De and north-south Dn deformation variables, and a complete three-dimensional deformation field De, Dn, Dz is obtained;

[0010] Step S4: Using the monitoring computing device, subsequent analysis is performed on the three-dimensional deformation field data obtained in step S3 to delineate the settlement boundary and identify the separation grouting target area;

[0011] Step S5: Using the monitoring computing device, based on the obtained settlement boundary and grouting target area data, grouting management evaluation decisions and hierarchical early warning rules are generated.

[0012] The specific method of step S1 is:

[0013] S101: A Delaunay triangular network is constructed by using the Delaunay triangular network method:

[0014] If it is C-band Sentinel-1 data, set the time baseline threshold to be ≤100 days and the spatial baseline threshold to be ≤150 meters;

[0015] If it is X-band TerraSAR-X data, set the time baseline threshold to be ≤80 days and the spatial baseline threshold to be ≤80 meters;

[0016] S102: Through differential interference, phase unwrapping, atmospheric correction and deformation sequence calculation, the average surface deformation rate field and the cumulative deformation variable of the entire mining area within the monitoring period are obtained;

[0017] S103: Based on the deformation rate field obtained in S102, set the deformation rate threshold to be -20 mm per year, identify the significant settlement area with a deformation rate ≤-20 mm per year, and mark the area as a key area that needs further fine monitoring. The boundary of the key area is marked by a vector polygon;

[0018] S104: Calculate the coherence coefficient of each pixel, set the coherence coefficient threshold to be 0.3, and mark the area with a coherence coefficient <0.3 as an InSAR incoherent area as a basis for subsequent unmanned aerial vehicle LiDAR aerial survey range planning to avoid invalid coverage.

[0019] The specific method of the step S2 is:

[0020] S201: The unmanned aerial vehicle flies according to a planned route, and high-density laser point cloud data is acquired;

[0021] The original point cloud data is denoised and classified, and a digital elevation model DEM_T1 of a T1 period and a digital elevation model DEM_T2 of a T2 period are generated based on ground points;

[0022] S202: Difference calculation is performed on the generated DEM_T1 and DEM_T2, and a vertical direction settlement amount Dz of a key area and an incoherent area is calculated, and the calculation formula is:

[0023] Dz = DEM_T2-DEM_T1;

[0024] Wherein, the unit of Dz is millimeter, and Dz < 0 indicates that the point has settlement in the T1-T2 period.

[0025] The specific method of the step S3 is:

[0026] S301: Data preprocessing and spatial registration are performed, including:

[0027] Data alignment: the InSAR deformation field data and the three-dimensional deformation field Dz data obtained by LiDAR difference are unified to the same geographic coordinate system;

[0028] Resolution unification: the bilinear interpolation method is used to resample the InSAR deformation field data with lower spatial resolution to the same grid scale as the three-dimensional deformation field Dz data;

[0029] Mask extraction: the vector diagram of the key area and the incoherent area generated in the step S1 is used as a mask to extract the target area to be solved from the registered data;

[0030] S302: A joint adjustment model based on a physical mechanism is constructed, including:

[0031] An observation equation is established:

[0032] The radar line-of-sight deformation variable Dlos of the InSAR deformation field data is calculated, and the calculation formula is:

[0033] Dlos = [cosθ×sinα]×De + [cosθ×cosα]×Dn + sinθ×Dz;

[0034] Wherein, De is an east-west direction deformation variable; Dn is a north-south direction deformation variable; Dz is a vertical direction deformation variable; θ is a known radar incidence angle, and α is a known satellite flight azimuth angle;

[0035] Introducing LiDAR constraint, the Dz value obtained by LiDAR in S202 is substituted into the above equation for solving;

[0036] Constructing the objective function:

[0037] Introducing the weight coefficient based on the observation accuracy, the following joint adjustment objective function is constructed:

[0038] ;

[0039] Wherein, n is the total number of effective pixels in the target area;

[0040] , are the InSAR measured line-of-sight deformation and LiDAR measured vertical subsidence of the i-th pixel, respectively;

[0041] , are the model predicted values calculated according to the current De, Dn estimated values and geometric relationship of the i-th pixel, respectively;

[0042] w los , w z are the weight coefficients of InSAR and LiDAR observations, respectively;

[0043] Define the LiDAR vertical measurement accuracy as σ z ≤3mm, and the InSAR line-of-sight measurement accuracy as σ los ≈15mm, then the weight ratio is set as:

[0044] w z / w los =(σ los / σ z )²≈(15 / 3)²=25;

[0045] To balance the influence of both, take w z =5, w los =1;

[0046] S303: Model solving and spatial continuity optimization, including:

[0047] Using weighted least squares method to iteratively solve the objective function constructed in S302, the optimal east-west deformation De and north-south deformation Dn are independently solved for each pixel;

[0048] Spatial filtering post-processing:

[0049] Introducing adaptive median filtering to smooth the De, Dn of the preliminary solution, and dynamically adjusting the filter window size according to the coherence coefficient of the pixel;

[0050] The obtained De and Dn are integrated with the known Dz of LiDAR to generate a complete, high-density and high-precision three-dimensional deformation field (De, Dn, Dz) of the target area;

[0051] S304: accuracy verification of the fusion result and model feedback, including:

[0052] A plurality of GNSS reference stations are arranged in the target mining area, and the GNSS measured three-dimensional deformation values (De_GNSS, Dn_GNSS, Dz_GNSS) are compared with the values of the three-dimensional deformation field at the corresponding positions obtained by fusion calculation;

[0053] Error calculation and threshold judgment:

[0054] The error of each direction deformation variable is calculated based on the difference between the fusion value and the GNSS measured value, and is defined as:

[0055] The absolute value of the error of the horizontal direction (De, Dn) is less than or equal to 8mm;

[0056] The absolute value of the error of the vertical direction (Dz) is less than or equal to 5mm;

[0057] If the threshold is exceeded, it indicates that the adjustment model is locally ill-posed;

[0058] When the verification error is out of limit, the system automatically starts the feedback optimization mechanism, and the optimization strategies include:

[0059] Re-evaluate and adjust the weight coefficients (w los , w z ) in S302;

[0060] In the local area with large error, spatial smoothing constraint is used to re-calculate the adjustment.

[0061] The specific method of the step S4 is:

[0062] S401: based on the vertical settlement Dz, Kriging interpolation is used to generate a settlement contour map, and the contour line with a cumulative settlement of 10mm is extracted as a settlement boundary;

[0063] S402: based on the comprehensive analysis of the three-dimensional deformation field (De, Dn, Dz), the following off-layer grouting target areas are identified:

[0064] The vertical deformation gradient mutation area with a Dz gradient greater than 5mm / 10m, which corresponds to the potential position of off-layer development;

[0065] The horizontal deformation significant area with a De or Dn deformation rate greater than 10mm / year, which reflects the surface horizontal stress concentration;

[0066] The above two types of areas are delineated as key target areas for off-layer grouting.

[0067] The specific method of the step S5 is:

[0068] S501: Grouting priority division: considering the following factors, the identified grouting target area is divided into high, medium and low priority:

[0069] Deformation gradient size, the larger the gradient, the higher the priority;

[0070] Horizontal deformation rate, the larger the rate, the higher the priority;

[0071] The distance from the key infrastructure in the mining area is closer, the priority is higher;

[0072] S502: Hierarchical early warning and grouting scheme linkage: obtain the average expansion speed V of the settlement boundary in the latest period and the minimum distance L from the infrastructure, set three early warning standards, and directly associate with the grouting scheme, The standard setting rule is:

[0073] When the obtained data is V<1 meter / month and L>200 meters, blue warning is carried out: at this time, the risk is low, and preventive grouting is recommended for high-priority grouting target area;

[0074] When the obtained data is 1 meter / month≤V≤3 meters / month or 100 meters≤L≤200 meters, yellow warning is carried out: at this time, the risk is medium, and monitoring needs to be strengthened, and control grouting is implemented for medium and high-priority grouting target area;

[0075] When the obtained data is V>3 meters / month or L<100 meters, red warning is carried out: at this time, the risk is high, and emergency response is needed immediately, and emergency treatment grouting is started for all identified grouting target area, especially high-priority target area;

[0076] S503: The system automatically generates a comprehensive report containing monitoring results, grouting target area distribution, priority, and treatment suggestion, and pushes it through safety management platform, short message, and on-site sound and light alarm.

[0077] The present application has the beneficial effects relative to the prior art: in view of the deficiencies of the existing monitoring technology in terms of dimension integrity and governance engineering connection, the present application aims to break through the limitations of traditional technical route, by deeply integrating the wide coverage ability of satellite remote sensing and the high-precision vertical measurement advantage of unmanned aerial vehicle LiDAR, a new type of monitoring method capable of accurately solving three-dimensional deformation field is established, and the monitoring results are directly applied to the target area identification, scheme design and treatment effect evaluation of delamination grouting, forming a integrated technical solution of monitoring-evaluation-treatment, and providing reliable technical support for the accurate prevention and control of mining area subsidence disaster. BRIEF DESCRIPTION OF DRAWINGS

[0078] The present application will be further described below in conjunction with the drawings:

[0079] Figure 1 The step flow chart of the mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion of the application. DETAILED DESCRIPTION

[0080] As Figure 1 shown, the application provides a mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion, specifically relates to a mine goaf subsidence precise monitoring method based on InSAR and unmanned aerial vehicle LiDAR cooperation and a corresponding grouting treatment evaluation method, solves the problem of mine area horizontal deformation monitoring, directly applies the obtained high-precision three-dimensional deformation field to the engineering guidance of separation layer grouting, forms a closed-loop technical system from precise monitoring to targeted treatment, realizes the technical scheme from wide-area screening, targeted encryption, model fusion, three-dimensional calculation to grouting decision, and the core scheme adopted includes:

[0081] Three-dimensional deformation field calculation model: the application takes the high-precision vertical deformation variable (Dz) obtained by the unmanned aerial vehicle LiDAR as a known constraint condition, and performs joint adjustment with the radar line-of-sight deformation variable (Dlos) obtained by the InSAR, so as to calculate the eastward (De) and northward (Dn) deformation variables which are difficult to directly measure, and finally obtain the specific algorithm model and method of the complete ground three-dimensional deformation field (De, Dn, Dz);

[0082] Grouting decision mechanism based on three-dimensional deformation field: the application identifies the deformation gradient mutation area and the horizontal deformation significant area as potential target areas of separation layer development based on the high-precision and high-density three-dimensional deformation data obtained by calculation, and generates a graded early warning and a differentiated grouting treatment scheme according to the deformation rate, the boundary expansion speed and the distance from the key facilities.

[0083] The application also provides a monitoring calculation device for realizing the corresponding functions, including a calculation device (such as a server or a ground workstation) with a processor and a memory.

[0084] To achieve the above purpose, the application provides a mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion, specifically including the following steps:

[0085] Step S1: wide-area subsidence screening and key monitoring area delineation based on SBAS-InSAR:

[0086] The computing device obtains time-series SAR image data covering the target mining area, and processes the data using SBAS-InSAR technology to obtain the deformation rate field and cumulative deformation of the entire mining area, identify significant subsidence areas and incoherent areas, and generate an overlay vector map of the key areas and incoherent areas. The computing device obtains time-series synthetic aperture radar (SAR) image data covering the target mining area, and preferentially selects C-band Sentinel-1 data (VV polarization, spatial resolution 5m x 20m) or X-band TerraSAR-X data (spatial resolution 1m x 1m); the SAR data is preprocessed, including thermal noise removal, radiation calibration, and multi-view processing (4 views in azimuth and 2 views in range), and then processed using SBAS-InSAR technology, with the following specific steps:

[0087] S101: Interferometric pair combination: Delaunay triangulation is used to construct an interferometric pair combination network, and if C-band Sentinel-1 data is used, the time baseline threshold is set to ≤100 days and the spatial baseline threshold is set to ≤150 meters; if X-band TerraSAR-X data is used, the time baseline threshold is set to ≤80 days and the spatial baseline threshold is set to ≤80 meters, to ensure that the interferometric pair quality meets the deformation inversion requirements.

[0088] S102: Deformation inversion: Through differential interference, phase unwrapping, atmospheric correction (using ERA5 atmospheric data for tropospheric delay correction), and deformation sequence solving, the average surface deformation rate field and cumulative deformation of the entire mining area during the monitoring period are obtained, with the deformation rate being in millimeters per year and the cumulative deformation being in millimeters.

[0089] S103: Key area delineation: Based on the deformation rate field, a deformation rate threshold of -20 millimeters per year is set (negative sign indicating subsidence), and significant subsidence areas with a deformation rate ≤-20 millimeters per year are identified and delineated as key areas that require further detailed monitoring, with the key area boundary marked by a vector polygon.

[0090] S104: Incoherent area marking: Calculate the coherence coefficient of each pixel, set a coherence coefficient threshold of 0.3, and mark areas with a coherence coefficient <0.3 as InSAR incoherent areas, generating an overlay vector map of the key areas and incoherent areas to provide a basis for subsequent UAV LiDAR survey range planning and avoid invalid coverage.

[0091] Step S2: Targeted detailed monitoring based on UAV LiDAR

[0092] The monitoring computing device uses the superimposed vector diagram of the focus area and the incoherent area generated according to step S1 to automatically plan the aerial survey range (completely covering the focus area and the incoherent area, and extending 50 meters outside the boundary to ensure data continuity) and the flight path (the flight path interval is set according to the point cloud density requirement to ensure that the grid accuracy of the generated digital elevation model (DEM) is better than 5 cm) of the unmanned aerial vehicle LiDAR; at two time points (T1 and T2 are separated by 6-12 months, and are adjusted according to the mining area subsidence rate) of the initial subsidence T1 and the subsidence development T2, the unmanned aerial vehicle LiDAR aerial survey task is performed, and the DJI Matrice350RTK unmanned aerial vehicle is used to carry the RIEGL VUX-1LR LiDAR device, and the specific steps are as follows:

[0093] S201: Data acquisition and processing: the unmanned aerial vehicle LiDAR system flies according to the planned flight path to obtain high-density laser point cloud data (point cloud density ≥ 20 points per square meter); the original point cloud data is processed, including denoising (removing noise points and isolated points), classification (distinguishing ground points and non-ground points such as vegetation and buildings), and generating DEM (denoted as DEM_T1) at T1 period and DEM (denoted as DEM_T2) at T2 period based on the ground points, and the grid accuracy of the DEM is better than 5 cm.

[0094] S202: Difference calculation: the DEM_T1 and the DEM_T2 are subjected to difference calculation to obtain the vertical subsidence Dz of the focus area and the incoherent area, and the calculation formula is:

[0095] Dz = DEM_T2-DEM_T1;

[0096] Wherein Dz is in millimeters, and Dz < 0 indicates that the point has subsided in the T1-T2 period.

[0097] Step S3: joint adjustment based on LiDAR vertical constraint, three-dimensional deformation field calculation and verification:

[0098] The monitoring computing device fuses the InSAR deformation data (average deformation rate field, cumulative deformation) obtained in step S1 and the LiDAR vertical subsidence Dz obtained in step S2 to calculate the east-west (De) and north-south (Dn) deformation, thereby obtaining the complete three-dimensional deformation field (De, Dn, Dz). The key of this step is to construct a joint adjustment model with high-precision Dz of LiDAR as a strong constraint, which solves the underdetermined problem of one equation and three unknowns in the InSAR observation equation. The specific steps are as follows:

[0099] S301: Data preprocessing and spatial registration:

[0100] Data alignment: unify the InSAR deformation field (usually raster data in geographic coordinate system) and the Dz field obtained by LiDAR difference (usually higher resolution raster data) to the same geographic coordinate system (WGS84 UTM coordinate system).

[0101] Resolution unification: the bilinear interpolation method is used to resample the InSAR deformation field data with lower spatial resolution to the same grid scale as the LiDAR Dz data (such as 1 meter grid) (this is to ensure that each registered pixel contains Dlos and Dz observation values at the same time, and to lay the foundation for pixel-by-pixel solving).

[0102] Mask extraction: use the "key area and out-of-coherence area" vector map generated in step S1 as a mask to accurately extract the target area to be solved from the registered full field data, greatly reducing unnecessary calculation and improving solving efficiency.

[0103] S302: Construct a joint adjustment model based on physical mechanism:

[0104] Establish observation equation: InSAR radar line-of-sight deformation Dlos and the true three-dimensional deformation component of the ground (De, Dn, Dz; De: east-west direction deformation variable, unit: millimeter; Dn: north-south direction deformation variable, unit: millimeter; Dz: vertical direction deformation variable, unit: millimeter) have the following strict geometric relationship:

[0105] Dlos=[cosθ×sinα]×De+[cosθ×cosα]×Dn+sinθ×Dz;

[0106] Wherein, θ is the radar incidence angle (known quantity), and α is the satellite flight azimuth angle (known quantity). This equation originally contains three unknowns De, Dn and Dz, which cannot be directly solved.

[0107] Introduce LiDAR constraint: the high-precision and high-density Dz value obtained by LiDAR in S202 is substituted into the above equation as a known quantity. At this time, the unknown number in the equation is reduced to two (De, Dn), so that the equation changes from "underdetermined" to "determined" or "overdetermined", which is mathematically solvable.

[0108] Construct objective function: considering that both InSAR observation value (Dlos) and LiDAR observation value (Dz) exist errors, and the precision is different (the precision of LiDAR is significantly higher than that of InSAR), the invention introduces a weight coefficient based on the accuracy of the observation value, and constructs the following joint adjustment objective function:

[0109] ;

[0110] Wherein, n is the total number of effective pixels in the target area.

[0111] , respectively are the InSAR measured line-of-sight deformation and the LiDAR measured vertical subsidence of the i-th pixel;

[0112] , respectively are the model predicted values of the i-th pixel calculated according to the current De, Dn estimated values and the geometric relationship;

[0113] w los , w z respectively are the weight coefficients of the InSAR and LiDAR observation values;

[0114] The determination basis is the ratio of the reciprocal square of the observation accuracy (root mean square error) of both.

[0115] Specifically, the LiDAR vertical measurement accuracy σ z ≤3mm, the InSAR line-of-sight measurement accuracy σ los ≈15mm, so the weight ratio can be set as:

[0116] w z / w los =(σ los / σ z )²≈(15 / 3)²=25;

[0117] In the embodiment, in order to balance the influence of both, w z =5, w los =1 (this quantitative weighting method based on error theory is an important guarantee for the scientificity and accuracy of the model, which avoids subjective experience assignment).

[0118] S303: Model solving and spatial continuity optimization:

[0119] Least square solving: the weighted least square method is used to iteratively solve the above objective function, and the optimal east-west deformation De and north-south deformation Dn are independently calculated for each pixel.

[0120] Spatial filtering post-processing: since De and Dn are independently solved for each pixel, there may be noise in space. Therefore, the adaptive median filtering is introduced to smooth the De and Dn field calculated initially. The filtering window size can be dynamically adjusted according to the coherence coefficient of the pixel (such as using a small window to retain details in a region with a high coherence coefficient, and using a large window to suppress noise in a region with a low coherence coefficient), while removing noise and maximizing the preservation of real deformation spatial details.

[0121] Three-dimensional deformation field generation: integrate the obtained De, Dn field with the known LiDAR Dz field to generate a complete, high-density, and high-precision three-dimensional deformation field (De, Dn, Dz) of the target area.

[0122] S304: fusion result accuracy verification and model feedback:

[0123] GNSS reference station verification: 3-5 GNSS reference stations are arranged in the target mining area, providing true value data independent of InSAR and LiDAR. Compare the GNSS measured three-dimensional deformation values (De_GNSS, Dn_GNSS, Dz_GNSS) with the values of the three-dimensional deformation field calculated by fusion at the corresponding position.

[0124] Error calculation and threshold judgment: calculate the error of deformation in each direction (fusion value-GNSS measured value). The absolute value of the horizontal direction (De, Dn) error is required to be ≤8mm, and the absolute value of the vertical direction (Dz) error is required to be ≤5mm. If the threshold is exceeded, it indicates that the adjustment model may be locally ill-posed.

[0125] Model feedback optimization: when the verification error exceeds the limit, the system automatically starts the feedback optimization mechanism. The optimization strategies include:

[0126] Re-evaluate and adjust the weight coefficients (w los , w z ) in S302;

[0127] In the local area with large error, use stronger spatial smoothing constraint to re-adjust the calculation.

[0128] This closed-loop verification mechanism ensures the reliability and accuracy of the final output three-dimensional deformation field data.

[0129] Step S4: settlement boundary demarcation and delamination grouting target area identification:

[0130] Use the high-precision three-dimensional deformation field data obtained in step S3 to perform subsequent analysis using the monitoring calculation device. The specific steps are as follows:

[0131] S401: Precise demarcation of settlement boundary: based on the vertical settlement Dz, use Kriging interpolation to generate a settlement contour map, and extract the contour line with a cumulative settlement of 10mm as the settlement boundary.

[0132] S402: Delamination grouting target area identification: based on the comprehensive analysis of the three-dimensional deformation field (De, Dn, Dz), identify two key areas:

[0133] Area One: Vertical deformation gradient mutation area (Dz gradient greater than 5mm / 10m), which often corresponds to the potential location of delamination development;

[0134] Region two: significant horizontal deformation area (De or Dn deformation rate greater than 10 mm / year), which reflects the surface horizontal stress concentration, is the key area of grouting management. The above two types of areas are delineated as the key target area of separation grouting.

[0135] Step S5: Grouting decision and hierarchical early warning:

[0136] Using the monitoring computing device, based on the settlement boundary and grouting target area, a comprehensive management decision is generated, including:

[0137] S501: Grouting priority division: considering the following factors, the identified grouting target area is prioritized (high, medium, low):

[0138] a) Deformation gradient size (the larger the gradient, the higher the priority);

[0139] b) Horizontal deformation rate (the larger the rate, the higher the priority);

[0140] c) Distance to key infrastructure in the mining area (such as buildings, roads, pipelines) (the closer the distance, the higher the priority).

[0141] S502: Hierarchical early warning and grouting scheme linkage: obtain the average expansion speed (V) of the settlement boundary in the last 6 months and the minimum distance (L) to the infrastructure, set three early warning standards, and directly associate with the grouting scheme:

[0142] Blue warning (V <1 meter / month and L > 200 meters): low risk, recommend preventive grouting for high-priority grouting target areas;

[0143] Yellow warning (1 meter / month ≤ V ≤ 3 meters / month or 100 meters ≤ L ≤ 200 meters): moderate risk, need to strengthen monitoring, and implement control grouting for high-priority grouting target areas;

[0144] Red alert (V > 3 meters / month or L < 100 meters): high risk, need immediate emergency response, and start emergency grouting for all identified grouting target areas, especially high-priority target areas.

[0145] S503: Output results and information push: the system automatically generates a comprehensive report containing monitoring results, grouting target area distribution, priority and management recommendations, and pushes it through multiple channels such as safety management platform, SMS, on-site sound and light alarm, etc.

[0146] The application breaks through the limitation of single monitoring technology, and realizes true three-dimensional deformation monitoring: the traditional satellite remote sensing technology cannot directly obtain the horizontal deformation component, and ground measurement is difficult to realize full coverage, the application combines the wide area coverage advantage of satellite remote sensing and the high-precision vertical measurement capability of unmanned aerial vehicle LiDAR through multi-source data fusion, effectively solves the underdetermined problem of traditional radar remote sensing one equation, three unknowns, realizes the simultaneous accurate calculation of east-west, north-south and vertical deformation of the mining area, and provides complete stereoscopic data support for mining area deformation analysis; based on the quantitative analysis of three-dimensional deformation field, the deformation gradient mutation area and horizontal stress concentration area can be accurately identified, these areas are the key parts of separation development, compared with the traditional experience-dependent grouting target area determination method, the target area identification provided by the application is more objective and accurate, effectively avoids the blindness of grouting engineering, improves the utilization efficiency of treatment resources, and provides a reliable decision basis for grouting engineering design; at the same time, by directly associating multi-period monitoring data with grouting decision, a response mechanism of settlement development trend and grouting emergency level is established, realizing the important change from passive monitoring to active intervention, when the settlement boundary accelerates expansion or deformation intensifies, the system can automatically improve the warning level and trigger the corresponding grouting emergency plan, significantly improving the initiative and response efficiency of mine safety management.

[0147] The application can objectively evaluate the treatment effect of grouting engineering by comparing the multi-period three-dimensional deformation field before and after grouting, monitor key indicators such as settlement rate change and deformation gradient improvement, this closed-loop verification mechanism not only provides a scientific evaluation method for grouting effect, but also feeds back the optimization of subsequent grouting parameter design, forms a continuous improvement technical cycle; in addition, the application realizes systematic and automatic processing of the whole process from data acquisition, processing, analysis to warning pushing and treatment suggestion generation, reduces the subjective error caused by manual intervention. The method is especially suitable for long-term monitoring demand in complex mining environment, significantly improves the work efficiency and timeliness of data, and provides stable and reliable technical support for mine safety production.

[0148] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the application, and not to limit them; although the application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the application.

Claims

1. A mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion, characterized in that: The monitoring evaluation step comprises the following steps: Step S1: using a monitoring computing device to obtain time-series SAR image data covering the target mining area, preferably C-band Sentinel-1 data or X-band TerraSAR-X data, performing thermal noise removal, radiation scaling, multi-view processing on the SAR image data, and then using SBAS-InSAR technology for processing to obtain InSAR deformation field data of the entire mining area, including deformation rate field and cumulative deformation, identifying significant subsidence areas and incoherent areas, and generating superimposed vector maps of key areas and incoherent areas; Step S2: using a monitoring computing device to automatically plan the aerial survey range and route of the unmanned aerial vehicle LiDAR according to the superimposed vector maps of the key areas and incoherent areas generated in step S1, and using the unmanned aerial vehicle to carry a laser radar device to perform unmanned aerial vehicle LiDAR aerial survey tasks at two time points, i.e., the initial subsidence period T1 and the subsidence development period T2, to obtain the vertical direction subsidence Dz of the key areas and incoherent areas; Step S3: using a monitoring computing device to fuse the InSAR deformation field data obtained in step S1 with the vertical direction subsidence Dz obtained in step S2 to calculate the east-west direction De and north-south direction Dn deformation, and obtain the complete three-dimensional deformation field De, Dn, Dz; Step S4: using a monitoring computing device to perform subsequent analysis on the three-dimensional deformation field data obtained in step S3 to delineate the subsidence boundary and identify the separation grouting target area; Step S5: using a monitoring computing device to generate grouting management evaluation decisions and hierarchical early warning rules based on the obtained subsidence boundary and grouting target area data.

2. The mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion according to claim 1, characterized in that: The specific method of step S1 is as follows: S101: using the Delaunay triangulation method to construct an interference pair combination network: If it is C-band Sentinel-1 data, set the time baseline threshold to be ≤100 days and the spatial baseline threshold to be ≤150 meters; If it is X-band TerraSAR-X data, set the time baseline threshold to be ≤80 days and the spatial baseline threshold to be ≤80 meters; S102: through difference interference, phase unwrapping, atmospheric correction, and deformation sequence calculation, obtain the average surface deformation rate field and cumulative deformation of the entire mining area within the monitoring period; S103: based on the deformation rate field obtained in S102, set the deformation rate threshold to be -20 mm per year, identify the significant subsidence area with a deformation rate ≤-20 mm per year, and mark this area as a key area that needs further fine monitoring, and the key area boundary is marked by a vector polygon; S104: calculate the coherence coefficient of each pixel, set the coherence coefficient threshold to be 0.3, and mark the area with a coherence coefficient <0.3 as an InSAR incoherent area as a basis for subsequent unmanned aerial vehicle LiDAR aerial survey range planning to avoid invalid coverage.

3. The mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion according to claim 2, characterized in that: The specific method of step S2 is as follows: S201: the unmanned aerial vehicle flies according to the planned route to obtain high-density laser point cloud data; perform denoising and classification processing on the original point cloud data, and generate a digital elevation model DEM_T1 for the T1 period and a digital elevation model DEM_T2 for the T2 period based on the ground points; S202: Difference calculation is performed on the generated DEM_T1 and DEM_T2 to obtain the vertical direction subsidence Dz of the key area and the incoherent area, and the calculation formula is as follows: Dz = DEM_T2 - DEM_T1; Wherein: the unit of Dz is millimeter, and Dz < 0 indicates that the point has subsidence in the period of T1-T2.

4. The mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion according to claim 3, characterized in that: The specific method of the step S3 is as follows: S301: Data preprocessing and spatial registration are performed, including: Data alignment: the InSAR deformation field data and the three-dimensional deformation field Dz data obtained by LiDAR difference are unified to the same geographic coordinate system; Resolution unification: the bilinear interpolation method is adopted to resample the InSAR deformation field data with lower spatial resolution to the same grid scale as the three-dimensional deformation field Dz data; Mask extraction: the vector diagram of the key area and the incoherent area generated in the step S1 is used as a mask to extract the target area to be solved from the registered data; S302: A joint adjustment model based on physical mechanism is constructed, including: An observation equation is established: The radar line-of-sight deformation Dlos of the InSAR deformation field data is calculated, and the calculation formula is as follows: Dlos = [cosθ×sinα]×De + [cosθ×cosα]×Dn + sinθ×Dz; Wherein, De is the east-west direction deformation variable; Dn is the south-north direction deformation variable; Dz is the vertical direction deformation variable; θ is the known radar incidence angle, and α is the known satellite flight azimuth angle; The Dz value obtained by LiDAR in S202 is introduced into the above equation for calculation; A target function is constructed: The weight coefficient based on the observation value accuracy is introduced to construct the following joint adjustment target function: ； Wherein, n is the total number of effective pixels in the target area; , respectively represent the InSAR measured line-of-sight deformation and the LiDAR measured vertical settlement of the i-th pixel. , are the model predicted values of the i-th pixel according to the current De, Dn estimates and the geometric relations, respectively; w los , w z are the weight coefficients of InSAR and LiDAR observations, respectively; The LiDAR vertical measurement accuracy is defined as σ z ≤ 3 mm, and the InSAR line-of-sight measurement accuracy is σ los ≈ 15 mm, the weight ratio is set as: w z / w los =(σ los / σ z )²≈(15 / 3)²=25; To balance both effects, take w z = 5, w los = 1; S303: Model solving and spatial continuity optimization are performed, including: The weighted least square method is adopted to iteratively solve the target function constructed in S302 to independently solve the optimal east-west direction deformation variable De and the south-north direction deformation variable Dn for each pixel; Spatial filtering post-processing: The adaptive median filtering is introduced to smooth the De and Dn solved initially, and the filter window size is dynamically adjusted according to the coherence coefficient of the pixel; The solved De and Dn are integrated with the known Dz of LiDAR to generate a complete, high-density and high-precision three-dimensional deformation field (De, Dn, Dz) of the target area; S304: Fusion result accuracy verification and model feedback, including: A plurality of GNSS reference stations are arranged in the target mining area, and the three-dimensional deformation values (De_GNSS, Dn_GNSS, Dz_GNSS) measured by GNSS are compared with the values of the three-dimensional deformation field at the corresponding positions solved by fusion; Error calculation and threshold judgment: The errors of the deformation variables in each direction are calculated based on the subtraction of the fusion values and the GNSS measured values, and the definitions are as follows: The absolute value of the horizontal direction (De, Dn) error is less than or equal to 8 millimeters; The absolute value of the vertical direction (Dz) error is less than or equal to 5 millimeters; If the threshold is exceeded, it indicates that the adjustment model is locally ill-posed; When the verification error is out of limit, the feedback optimization mechanism is automatically started, and the optimization strategies include: re-evaluating and adjusting the weight coefficients (w los , w z ) in S302; In the local area with large error, the spatial smoothing constraint is used to re-adjust the error.

5. The mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion according to claim 4, characterized in that: The specific method of step S4 is: S401: Based on the vertical settlement Dz, the Kriging interpolation is used to generate the settlement contour map, and the contour line with a cumulative settlement of 10 mm is extracted as the settlement boundary; S402: Based on the comprehensive analysis of the three-dimensional deformation field (De, Dn, Dz), the following off-layer grouting target areas are identified: The vertical deformation gradient mutation area with a gradient greater than 5 mm / 10 m, which corresponds to the potential position of off-layer development; The horizontal deformation area with a deformation rate greater than 10 mm / year, which reflects the surface horizontal stress concentration; The above two types of areas are delineated as key target areas for off-layer grouting.

6. The mine area subsidence monitoring and grouting treatment evaluation method based on multi-source remote sensing fusion according to claim 5, characterized in that: The specific method of step S5 is: S501: Grouting priority classification: considering the following factors, the identified grouting target areas are classified into high, medium and low priority: Deformation gradient, the larger the gradient, the higher the priority; Horizontal deformation rate, the larger the rate, the higher the priority; The distance to the key infrastructure in the mining area, the closer the distance, the higher the priority; S502: Hierarchical early warning and grouting scheme linkage: obtain the average expansion speed V of the settlement boundary in the latest period and the minimum distance L to the infrastructure, set three early warning standards, and directly associate with the grouting scheme, the standard setting rules are: When the obtained data is V<1 meter / month and L>200 meters, blue warning is given: at this time, the risk is low, and preventive grouting is recommended for high-priority grouting target areas; When the obtained data is 1 meter / month≤V≤3 meters / month or 100 meters≤L≤200 meters, yellow warning is given: at this time, the risk is medium, monitoring needs to be strengthened, and control grouting is implemented for medium and high-priority grouting target areas; When the obtained data is V>3 meters / month or L<100 meters, red warning is given: at this time, the risk is high, emergency response is needed immediately, and emergency treatment grouting is started for all identified grouting target areas, especially high-priority target areas; S503: The system automatically generates a comprehensive report containing monitoring results, grouting target area distribution, priority, and treatment recommendations, and pushes it through multiple channels such as safety management platform, SMS, and on-site sound and light alarm.

Images