Precise monitoring method for ground surface settlement of goaf grouting

Abstract

The invention discloses a goaf grouting ground surface settlement precise monitoring method, which uses a ground surface deformation precise monitoring method fusing GNSS, InSAR, D-InSAR and SBAS-InSAR technologies, generates an interference image set on the basis of InSAR and D-InSAR, selects a stable coherent target point, solves the phase information of the coherent target point, and obtains the goaf grouting ground surface settlement precise monitoring result. Therefore, information such as ground subsidence rate and time sequence of the coherent target point is obtained. According to the method, satellite remote sensing, a low-altitude unmanned aerial vehicle and a laser LiDAR are fused to assist a ground precision level construction site residual subsidence deformation air-space-ground integrated monitoring technology, efficient acquisition and refined monitoring of construction site ground surface deformation information are realized, a construction site refined model based on the unmanned aerial vehicle and the LiDAR is constructed, and the construction site residual subsidence deformation monitoring method based on the InSAR and the D-InSAR is established based on the InSAR and the D-InSAR. Generating an interference image set; the method is high in monitoring precision, wide in monitoring range, small in weather influence, convenient and easy in monitoring implementation, relatively low in cost, high in safety, fast in data updating and rich in data volume, and can monitor the continuous deformation process of a ground target on a time sequence.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of mining area monitoring, in particular to a goaf grouting surface subsidence precision monitoring method. BACKGROUND

[0002] Due to underground mining, the overburden strata collapse, fracture or bending, although the natural compaction effect for a long time, the mined-out area, fissure, separation layer and the like formed by mining will still exist for a long time, and the surface of the old goaf will still have a certain degree of residual subsidence deformation, which will cause damage and deformation to the newly built structures and equipment foundations. Especially for the goaf formed by strip mining, room-and-pillar mining, high-water filling mining, fault pillar mining and the like, the coal pillar is prone to instability, resulting in a large residual subsidence deformation of the surface. Goaf grouting and filling is an effective technical measure. The main mechanism of goaf grouting and filling is to occupy most of the old goaf stratum residual pore fissure, cavity and separation layer space by grouting body, prevent the residual coal pillar from instability, and control the residual deformation of the old goaf. However, goaf grouting cannot guarantee that the overburden strata will not settle again. Influenced by factors such as coal seam occurrence, mining method, stratum lithology, filling process, filling material and the like, the surface of the goaf after grouting and filling may still have relatively obvious movement and deformation. Moreover, the buildings above the goaf are now developing towards large-scale, with higher height and larger plane span, and the like, and such building load has a greater influence depth, and is more likely to cause the activation of the goaf stratum system. At present, there is less research on the surface movement law under the condition of goaf grouting and filling, and there is no systematic laboratory analysis and theoretical research, and there is a lack of theoretical basis and calculation model for reference.

[0003] Therefore, the person skilled in the art provides a goaf grouting surface subsidence precision monitoring method to solve the problems in the background art. SUMMARY

[0004] To solve the above technical problems, the present application provides a goaf grouting surface subsidence precision monitoring method, characterized in that it comprises the following steps:

[0005] I. Establish a ground subsidence monitoring control network; research on grouting and filling of goaf surface precision leveling observation technology Establish a ground subsidence monitoring control system, lay out a high-precision monitoring control network according to the second-order leveling measurement accuracy, monitor and research the surface residual subsidence under the condition of goaf grouting, adopt overall observation of measuring points combined with key point review, and construct a point-surface combined subsidence monitoring network to ensure the field observation accuracy and improve the data acquisition efficiency; establish an air-space monitoring network system; it integrates satellite remote sensing, low-altitude unmanned aerial vehicle and laser Lidar technology to collect data.

[0006] II. Observation data processing and analysis; Through InSAR technology, laser radar point cloud, optical remote sensing image processing technology and other technical means, the ground subsidence and land use and cover change of the study area are obtained, the adjustment model is optimized, the field data processing is carried out, and the surface movement deformation characteristics of the grouting influence area are analyzed, which provides measured basis for the study of surface residual deformation law of grouting filling goaf.

[0007] III. Similar material simulation experiment: using similar material simulation experiment, numerical simulation calculation, mechanical analysis and other methods, the movement law of overburden strata in grouting filling goaf is studied, and the movement mechanism of overburden strata in grouting filling goaf is revealed

[0008] IV. Experimental analysis: experimental test of physical and mechanical properties of soil, rock, filling body and filling body and surrounding rock composite, obtain rock mechanics characteristics and parameters of goaf under different confining pressures and unloading rates, analyze and evaluate the mechanical properties of overburden strata before and after grouting filling in goaf; take the caving zone rock and the filling body in the goaf during construction, carry out hierarchical loading creep test, obtain the creep curve and long-term strength change law, establish the corresponding constitutive model and carry out parameter identification;

[0009] V. Through the processing and analysis of observation data and simulation experiment data, the coupling relationship between residual deformation of goaf and compression deformation of building foundation is analyzed.

[0010] Preferably, in step one, the space-air monitoring network system is a space-air integrated monitoring technology system based on Beidou, which integrates GNSS and InSAR technology for precise and accurate monitoring of ground deformation. The system combines the advantages of GNSS (Beidou high-precision positioning) and InSAR technology to achieve precise and accurate real-time monitoring of ground deformation in the grouting engineering area, and to obtain the three-dimensional deformation field of the study area, providing data support for analyzing the surface deformation law of the grouting engineering area under the influence of multiple factors. The system uses unmanned aerial vehicle + LiDAR regional fine surface model, uses unmanned aerial vehicle carrying high-precision LiDAR to construct continuous high-precision DEM and DSM model, directly describes and expresses the terrain surface, and obtains the overall deformation result and deformation trend of the monitoring area. The system uses remote sensing technology to monitor the surface cover change and historical evolution law of ground subsidence, based on multi-temporal high-resolution remote sensing data, accurately extracts the historical evolution data of land, vegetation and other surface covers in the study area, and obtains the evolution law. The system uses SBAS-InSAR and other data processing technologies to obtain the average deformation rate information and deformation history information of ground targets in the study area.

[0011] Preferably, in step two, the method of combining D-InSAR and coherent target time series analysis is used to process the selected satellite radar data. The InSAR data processing includes data collection, radar data registration, D-InSAR data processing and effective data selection, data quality evaluation, etc. Specifically, the following steps are included:

[0012] 1) Data collection: In the work, the necessary InSAR data, DEM data and basic geographic information data for mapping analysis are collected;

[0013] 2) Radar data registration: Obtain the precise orbit data of each period of data, use orthogonal correlation technology to realize the precise registration between each main and auxiliary image, and the registration accuracy of the whole scene data is required to be higher than 0.2 pixels. Register the DEM data with the RadarSAT-2 main image as the reference, remove the terrain phase in subsequent interferometric measurement through DEM data information and radar imaging geometric information. Obtain the SRTM DEM data covering the test area, and according to the SAR image imaging geometry and other information, perform coordinate system conversion, i.e. geographic coding, on the external DEM, establish the conversion lookup table of the coordinate systems between the SAR image and the external DEM, and through the texture matching of the DEM simulated SAR image and the SAR image, the matching accuracy is higher than 0.2 pixels, and the lookup table is corrected. Through the lookup table, the DEM can be sampled to the SAR radar imaging coordinate system, and then combined with the satellite orbit information to simulate the correct terrain phase, and complete the difference processing of the interferometric phase. Similarly, in the later result output, the deformation information in the SAR imaging coordinate system can be sampled to the DEM geographic coordinate system, and the analysis and mapping work of the later results can be completed;

[0014] 3) Selection of interferometric image pairs and D-InSAR processing: Add precise orbit data to each period of data to reduce the noise and terrain error caused by orbit error and other influences;

[0015] 4) Coherent target time series analysis: Time series analysis is performed on the interferometric phase of the coherent point target, and according to the time and space characteristics of each phase component, the atmospheric fluctuation, DEM error and noise are estimated and separated from the differential interferometric phase one by one, and finally the deformation rate of each coherent point target is obtained. The accuracy of the annual deformation rate obtained can reach millimeter level;

[0016] Using the differential interferometric phase obtained in the early stage, the interferometric phase of the coherent point target is extracted, the spatial phase unwrapping is performed on the point target phase, the Delaunay triangular network is constructed according to the positions of the points, the residual points in the triangular network are judged, and the minimum cost flow algorithm is used to connect the positive and negative residual point pairs to obtain the unwrapping result of each point;

[0017] For the differential interference phase sequence of the coherent point target, the two-dimensional periodic function is used for parameter inversion to estimate the adjacent point moving speed and height error correction. Through multiple error correction iterations, the correlation coefficient of the two-dimensional periodic function is 0.8-0.85, and the linear deformation rate is accurately extracted;

[0018] 5) Data comprehensive analysis and mapping;

[0019] 6) Quality control of data processing: InSAR processing is very strong in process, and 100% human-computer interaction check and mutual check between operators are carried out at each step.

[0020] Preferably, in step 1), data collection, the specific steps are as follows:

[0021] (1) The radar data used in the InSAR ground subsidence information extraction work is taken from the RadarSAT-2 satellite, and the precise orbit parameters of each period of data are obtained;

[0022] (2) The DEM data involved in the InSAR data processing work is collected, and the DEM data obtained by SRTM is used, with a horizontal resolution of 30 meters, which meets the requirements of terrain compensation in InSAR processing;

[0023] (3) Other optical remote sensing image data, basic geographic information data and the like for assisting ground subsidence status analysis, result mapping and the like.

[0024] Preferably, in step 3), the selection of the interference image pair and the D-InSAR processing, the specific steps are as follows:

[0025] (1) Selection of interference image pairs: according to the baseline distribution characteristics of the obtained SAR data sequence, combined with the basic characteristics of regional surface deformation, each data and the adjacent 5 data form an interference image pair sequence, 180 interference image pairs are formed from 39 data, according to the 24-day data time sampling, the time baseline is between 24-120 days, due to the missing of the middle part of the data, the time baseline of part of the interference pairs is longer, and all are less than 240 days;

[0026] (2) Generation of interference fringe sequence: interference processing is performed on the interference image pairs that meet the baseline requirements to generate an interference fringe sequence, and adaptive filtering processing is performed on the interference fringe according to the characteristics of the noise phase;

[0027] (3) Differential interference processing: precise orbit parameters are introduced for baseline estimation, flat ground is removed according to the geometric imaging model, DEM sampled to the radar coordinate system is used, and the terrain phase is simulated through the geometric imaging model to obtain the differential interference graph containing deformation phase, atmospheric phase and noise phase;

[0028] (4) Interference fringe screening: screening the differential interference fringe, excluding the interference graph with obvious strong convection influence and the interference graph with obvious residual interference phase fringe.

[0029] Preferably, in step 5), the specific steps of data comprehensive analysis and mapping are as follows:

[0030] (1) The point target deformation rate information is sampled to the commonly used coordinate system through the coordinate lookup table created during DEM registration;

[0031] (2) The InSAR ground subsidence information is analyzed by comprehensively applying image processing and GIS technology, and the InSAR monitoring ground subsidence map is prepared.

[0032] Preferably, in step 6), the quality control of data processing specifically includes the following contents:

[0033] (1) Raw data quality: The overall inspection is carried out on the data type, coverage range, phase, beam mode, etc. of each radar data obtained for the goaf grouting and filling ground surface movement precision monitoring and residual deformation law research, the data is read into the professional processing software, the reading is correct, there is no bad point in the map, and the phase information is not lost;

[0034] (2) Image registration (radar data and DEM data): The radar data is registered by selecting the same main image, the registration accuracy is higher than 0.2 pixels in distance and direction, and the registration accuracy of DEM and main image reaches 0.2 pixels. If the registration accuracy requirement is not met, the effective phase cannot be formed in the subsequent differential interference processing;

[0035] (3) Differential interference processing: The differential interference processing includes interference of two SAR data, terrain phase removal, flat phase removal, and each differential interference phase graph formed is checked to ensure the accuracy of the differential coherence;

[0036] (4) Residual error correction: The residual terrain phase removal and orbit error removal are carried out on the differential phase graph formed after the differential interference processing, and the interference phase with obviously significant atmospheric error is excluded, and the next step calculation is not carried out;

[0037] (5) Coherent target point selection: The accuracy of the coherent point target position is checked, and whether the coherent point target density meets the subsequent time analysis requirement is checked. The point density in the urban area is required to be more than 100 per square kilometer;

[0038] (6) Coherent target time series analysis: The coherent point target time series analysis is carried out, and the ground subsidence rate information is extracted. The phase unwrapping file, residual error file and deformation file generated by the time series analysis are checked to ensure that each layer is unwrapped correctly, there is no discontinuous phase, and the residual error file has no blocking phenomenon;

[0039] (7) Geocoding: Geocoding deformation information, checking the correctness of the file geographic location;

[0040] (8) Format conversion: Convert deformation information into shp data format for analysis and mapping, check the correctness of shp geometry and attributes.

[0041] Preferably, in step two, the laser radar point cloud data is preprocessed, and the profile section method is used in data processing to check data quality. The specific steps are as follows:

[0042] 1) Render the point cloud into different colors according to the elevation value, that is, point cloud coloring. Point cloud coloring is to assign actual texture color to the original LiDAR point cloud data without RGB color information. The ground object information in the survey area can be intuitively displayed in the point cloud view. In the later stage, RGB color information can be used to assist point cloud classification. There are mainly two forms to realize point cloud coloring:

[0043] (1) Color according to the original image taken in the survey area, camera calibration file, POS file after aerial triangulation, and coloring range;

[0044] (2) Color according to the orthophoto generated after aerial triangulation or fast mosaic. The data processing uses the orthophoto-based method to realize point cloud coloring;

[0045] 2) Use the profile function to make a profile in the overlapping area of each flight strip in the main interface. Observe the profile view to see if there is obvious layering. If it is found that the point cloud is layered, it needs to be reprocessed for point cloud calculation and flight strip adjustment. In addition, the height difference of point cloud between flight strips can also be checked through the quality check function in the intelligent laser module of the unmanned aerial vehicle manager, which can also reflect whether the point cloud is layered from the side.

[0046] 3) Coordinate transformation. The point cloud coordinate system needed is CGCS2000 national geodetic coordinate system, Gauss 3-degree band projection, central meridian 117°, 1985 national height datum. After point cloud adjustment and coloring, coordinate transformation is performed. It mainly involves parameter calculation and coordinate transformation. The residuals dN, dE and dU of the seven parameters to be solved must be controlled within the limit. The detailed steps are as follows:

[0047] (1) The initial point cloud is in the WGS84 system. Now we need to get the point cloud in the CGCS2000 national geodetic coordinate system. We need to perform parameter conversion between the two coordinate systems. Use three parameter points measured in the field. At the same time, know the three pairs of coordinates of the three points in the WGS84 coordinate system and the CGCS2000 national geodetic coordinate system. Perform parameter calculation to get the seven parameters of the two coordinate systems. Check if the residual information dN, dE and dU meets the requirements;

[0048] (2) Coordinate transformation parameter configuration of the survey area is performed, the Bursa seven parameters obtained by parameter solving are added, the conversion parameter configuration is completed, and the point cloud is converted from the source coordinate system (WGS84 coordinate system) to the target coordinate system (CGCS2000 national geodetic coordinate system).

[0049] 4) Standard point cloud data output After the above preprocessing operation of the unmanned aerial LiDAR point cloud data, the point cloud data for production can be generated, and the obtained point cloud results are generally stored in the form of “*.las” format. The ground point cloud data and image data can be used to generate three-dimensional models, DEM, DSM and DOM digital products according to needs.

[0050] Preferably, after the point cloud data preprocessing, the obtained point cloud data needs to be further processed according to actual needs, and the detailed steps are as follows:

[0051] 1) Point cloud denoising

[0052] The obtained mine area standard format las point cloud can be visualized to obtain a large number of noise points, and there are noise point groups above and below the ground point cloud. In order to ensure the accuracy of the DEM data product generated by using the point cloud data, it is necessary to remove the noise point group, and the point cloud denoising function in the point cloud module of the unmanned aerial vehicle manager is used to denoise the point cloud data of the whole project area;

[0053] 2) Point cloud filtering and classification

[0054] The filtering of point cloud refers to the process of removing non-ground points (vegetation, buildings, bridges, power facilities, etc.) in point cloud data in order to extract digital elevation model from point cloud data. This process is called filtering of laser radar data.

[0055] The classification of point cloud refers to the process of classifying non-ground points, which is used for feature extraction and three-dimensional reconstruction of buildings.

[0056] 3) Digital product generation

[0057] After point cloud denoising and automatic classification of standard point cloud by using the point cloud module of the unmanned aerial vehicle manager, high-precision laser point cloud data for producing digital products is obtained. The data editing function in the point cloud module can be used to generate DSM from the denoised point cloud, and DEM can be further generated from the ground point cloud separated by automatic classification. For the area with poor automatic classification effect, manual classification can be performed through regional smoothing, online classification and offline classification. The DSM and DEM digital products generated by the point cloud obtained after denoising and classification processing. After point cloud denoising and automatic classification, DSM is generated by using the obtained ground point cloud.

[0058] 4) Data accuracy check mainly includes

[0059] Point cloud data density and point cloud data plane and elevation accuracy check, with unmanned aerial vehicle manager in the wisdom of laser module point density measurement can be measured on the point density of point cloud data, you can set the measurement box length value to measure the point density in different areas, selected three types of regional project area, respectively, on the point density measurement, three types of areas are building dense area, vegetation dense area, ground area, with accuracy check function can check the point cloud elevation accuracy, by importing the check point file, and by point cloud solution obtained by point cloud data, get the elevation difference between point cloud data and check point coordinates, realize the accuracy check of point cloud data.

[0060] Preferably: analyze the coupling relationship between residual deformation of goaf and compression deformation of building foundation, the detailed steps are as follows:

[0061] 1) Residual ground surface deformation calculation method under goaf grouting condition

[0062] (1) Establish the theoretical model of probability integral method

[0063] Goaf grouting filling can slow down the residual deformation of goaf ground surface, and the ground surface movement form is relatively gentle residual subsidence, which presents a basin-like shape with large middle and small edge above the goaf, and the subsidence basin form still conforms to the normal curve distribution law, so the probability integral method mathematical model can be used for residual ground surface deformation calculation. The theoretical model of probability integral method is established, and through mathematical derivation, the expression of the subsidence value We(x,z) of (x,z) point caused by unit mining can be obtained.

[0064]

[0065] In the formula:

[0066] W e (x,z)——subsidence value caused by unit mining, unit mm;

[0067] r z ——main influence radius, unit m,

[0068]

[0069] For the ground surface, z is equal to the mining depth, and A is a constant representing the medium unit size, so that r z is a constant, and let r z be a constant r, then the expression of the subsidence basin of the ground surface unit is:

[0070]

[0071] (2) Residual ground surface deformation calculation of goaf grouting filling

[0072] The subsidence of any point (x, y) on the surface is:

[0073]

[0074] wherein:

[0075] l A 、l B Length and width of the calculation area of the goaf, unit: m;

[0076] M d Mining thickness, unit: m;

[0077] q c Residual surface subsidence coefficient of the goaf grouting and filling;

[0078] α——Coal seam inclination;

[0079] r——Main influence radius, unit: m.

[0080] The inclination, curvature, horizontal movement and horizontal deformation of any point on the surface in the given φ direction are:

[0081]

[0082] (3) Analysis of residual surface subsidence calculation parameters

[0083] The residual surface deformation of the goaf grouting and filling can be calculated by the residual subsidence coefficient calculated by the goaf filling rate, and the calculation formula is as follows:

[0084] (a) Residual surface subsidence coefficient qc of the goaf grouting and filling

[0085] The residual surface subsidence coefficient is mainly related to the subsidence coefficient in the surface movement duration and the duration after the end of mining. The calculation method of the residual surface movement and deformation subsidence coefficient is as follows:

[0086]

[0087] q——Subsidence coefficient in the surface movement duration;

[0088] q1——Residual surface subsidence coefficient under general mining conditions;

[0089] k——Adjustment coefficient, generally 0.5-1.0;

[0090] t——Time from the end of mining;

[0091] The surface movement and deformation under the condition of goaf grouting and filling can be calculated by the method of calculating residual subsidence coefficient by goaf filling rate, and the residual surface subsidence coefficient of goaf grouting and filling can be calculated by the following formula:

[0092]

[0093] In the formula: q c residual surface subsidence coefficient of goaf grouting and filling;

[0094] ρ c goaf filling rate;

[0095] ρ y compaction rate of goaf filling body.

[0096] The compaction rate of goaf filling body can be measured by compression test, and the calculation method can be calculated by formula (11), and generally takes about 0.1

[0097]

[0098] In the formula:

[0099] V - filling body volume filled in goaf, unit m 3 ;

[0100] V' - filling body volume after compression, unit m 3

[0101] (b) inflection point offset coefficient S / H

[0102] When the overburden type is hard, medium hard, and weak, the S / H value when the surface is basically stable is 0.31-0.43, 0.08-0.30, and 0-0.07, respectively.

[0103] The inflection point offset coefficient when calculating the residual surface subsidence of goaf can be analogized to the parameters when the surface is basically stable.

[0104] (c) main influence angle tangent tanβ

[0105] The angle tangent tanβ does not change significantly with the extension of time after the end of mining, and its value is referred to the tanβ of conventional subsidence deformation, and when the overburden type is hard, medium hard, and weak, the value is 1.2-1.91, 1.92-2.40, and 2.41-3.54, respectively.

[0106] (d) horizontal movement coefficient b

[0107] The horizontal movement coefficient b is mainly related to the properties of overburden, and has little relation with other factors. The horizontal movement coefficient b in the calculation of residual surface subsidence can refer to the value of conventional subsidence deformation. The value range of the horizontal movement coefficient is 0.2-0.4, and generally is 0.3.

[0108] (e) Mining influence propagation angle θ

[0109] The mining influence propagation angle θ is dependent on the coal seam inclination and the mining influence propagation coefficient K. The mining influence propagation coefficient K is mainly related to the overburden lithology, coal seam mining depth, etc. The grouting and filling of goaf has limited influence thereon. The mining influence propagation angle θ in the calculation of residual surface subsidence can refer to the value of conventional subsidence deformation θ. θ = 90°-kα, α is the coal seam inclination, and the value of k is 0.7-0.8, 0.6-0.7, and 0.5-0.6 respectively when the overburden types are hard, medium hard, and weak.

[0110] 2) Calculation model of ground soil compression settlement under building load

[0111] When the engineering construction is carried out above the goaf, the newly built building (structure) will be affected by the compression settlement deformation of the ground soil under the building load in addition to the residual deformation of the goaf. Therefore, the residual surface deformation above the goaf should include the superposition of the residual movement deformation of the goaf and the compression settlement of the soil under the load.

[0112] (1) Spatial stress analysis of ground soil under load

[0113] The influence range of the rectangular load (a x b) of the building on the underlying ground soil is approximately an ellipsoid. For the convenience of calculation, the maximum influence range of the soil can be regarded as a cuboid with the length la, the width lb, and the height Hz. la and lb are the length and width of the soil affected by the load, and Hz is the load influence depth,

[0114] Therefore, the length of la and lb is taken as the connecting line of the critical points at which the load influence depth is 0 along the center lines of the long and short sides of the ground surface. The maximum value of Hz is taken as the load influence depth. A point M is taken in the affected space range of the soil. A unit body is taken at the point M under the uniform load P of the rectangular foundation. The additional stresses in the vertical and horizontal directions at the unit body can be calculated according to formula 12.

[0115]

[0116] In the formula: σ fz(x,y,z) , σ fx(x,y,z) , σ fy(x,y,z) Additional stresses in the x, y, and z directions at the point (x, y, z);

[0117] σ z(x,y,z)σ x(x,y,z) σ y(x,y,z) σ (X,Y,Z)

[0118] (2) Layered Settlement Model of Foundation Soil

[0119] The layered summation method in soil mechanics is applied to calculate the settlement of foundation considering the compression characteristics of the soil and the layered characteristics of the natural soil. The following assumptions are made in this model:

[0120] (a) Each layer in the foundation is compressed vertically without lateral expansion.

[0121] (b) The thickness of each layer is not too large, usually 1-2m, and the stress at the top and bottom of each layer is not constant, so the average stress is used in the calculation.

[0122] (3) Calculate the deformation of the soil within the load influence range.

[0123] The settlement of a point (x, y, 0) on the surface of the foundation soil is the cumulative compression settlement of the soil elements along the straight line from this point to the depth of z = Hz. The settlement value of the element at point M (x, y, z) can be calculated by formula 13, and dS (X,Y,Z) The integral is performed on z = [0, Hz], and the settlement at point (x, y, 0) can be calculated by formula 14:

[0124]

[0125]

[0126] In the formula:

[0127] dS (X,Y,Z) The vertical settlement of the element soil at point (x, y, z);

[0128] E S(X,Y,Z) The compression modulus of the element soil at point (x, y, z);

[0129] S (x,y,0) The settlement at point (x, y) on the surface of the foundation soil (z = 0).

[0130] Since the foundation soil has the characteristic of natural soil layering, from the plane z = 0 to the plane z = Hz, it can be divided into m1 soil layers, so S (x,y,0) It can also be calculated by formula 14:

[0131]

[0132] In the formula:

[0133] — the average additional stress of the i-th soil layer under the point (x, y, 0) ;

[0134] — the average additional stress coefficient of the i-th soil layer under the point (x, y, 0) ;

[0135] Δz — the thickness of the i-th soil layer under the point (x, y, 0) ;

[0136] E s — the compression modulus of the i-th soil layer under the point (x, y, 0) ;

[0137] m1 — the number of soil layers in the load influence depth range;

[0138] p — the average pressure of the foundation, unit kg / m2.

[0139] The settlement of the building foundation under the load can be considered as the superposition of the compression settlement of the foundation soil body and the residual surface subsidence, i.e. 总 = S d + S w In the settlement calculation of the two, the origin of the selected coordinate system is different, the compression calculation coordinate system of the foundation soil body is converted to the residual deformation calculation coordinate system of the goaf, and then the cooperative relationship formula of the residual settlement of the goaf and the compression deformation of the building foundation can be calculated according to formula 15.

[0140]

[0141] The settlement and horizontal movement of the foundation soil body are mainly the cooperative settlement of the residual deformation of the goaf and the compression deformation of the foundation soil body, and then the inclination, horizontal movement, curvature and horizontal deformation of a given p direction on the surface of the foundation soil body can be represented by the following formula

[0142]

[0143] (4) Building foundation settlement deformation calculation under the condition of goaf grouting filling

[0144] The settlement of the building foundation under the condition of goaf grouting filling mainly includes the residual surface subsidence of the goaf and the compression settlement of the foundation soil body, and the specific calculation is shown in formula 17-20. In the calculation, the residual surface subsidence coefficient can be calculated according to formula 9, and the compression settlement of the foundation soil body is calculated by using the layer summation method, so that the settlement of 1 / 4 area of the foundation soil body in the cuboid space can be calculated. The settlement of a point on the surface of the foundation soil body is the cumulative compression settlement of each unit soil body on the straight line from the point to the load influence depth, which needs to be layered according to the type and thickness of the soil body. The layering thickness should not be too large (not more than 3m), and the settlement of any point on the surface of the soil body is calculated by using formula 14.

[0145] Technical effects and advantages of the present application:

[0146] Compared with the prior art, the present application has the beneficial effects that:

[0147] 1、The present application is an air-space-ground integrated monitoring technology for residual subsidence deformation of construction site assisted by satellite remote sensing, low-altitude unmanned aerial vehicle and laser LiDAR, which realizes efficient acquisition and fine monitoring of surface deformation information of the construction site, constructs a fine model of the construction site based on unmanned aerial vehicle+LiDAR, generates an interferogram set based on InSAR and D-InSAR, has high monitoring precision, wide monitoring range, small weather influence, convenient and easy monitoring implementation, relatively low cost, high safety, fast data update, rich data quantity, and can monitor the continuous deformation process of the ground target in the time sequence.

[0148] 2、The present application analyzes the local uplift and subsidence deformation characteristics of the ground surface under the condition of goaf grouting and filling, establishes a ground subsidence deformation calculation model based on the coupling of residual deformation of the goaf and compression deformation of the building foundation, and deduces a calculation formula of residual subsidence coefficient of the ground surface under grouting and filling, which is beneficial to the processing and analysis of collected data and makes the monitoring result more accurate. BRIEF DESCRIPTION OF DRAWINGS

[0149] Figure 1 It is a flow chart of the goaf grouting ground settlement precision monitoring method provided by the present application.

[0150] Figure 2 It is a spatial relationship diagram of the goaf calculation area and the ground surface in the goaf grouting ground settlement precision monitoring method provided by the present application.

[0151] Figure 3 It is a spatial influence range diagram of the foundation soil body in the goaf grouting ground settlement precision monitoring method provided by the present application.

[0152] Figure 4 It is a spatial relationship diagram of the length l a , the width l b and the load influence depth H z of the soil body affected by the load in the goaf grouting ground settlement precision monitoring method provided by the present application.

[0153] Figure 5 It is a plane coordinate conversion relationship diagram of the foundation soil body compression settlement and residual ground subsidence in the goaf grouting ground settlement precision monitoring method provided by the present application. Figure 6 It is a medium particle physical model diagram in the goaf grouting ground settlement precision monitoring method provided by the present application. Figure 7 It is a medium particle movement probability model diagram in the goaf grouting ground settlement precision monitoring method provided by the present application. Figure 8 is a surface subsidence probability distribution map in a goaf grouting surface subsidence precision monitoring method provided by the embodiment of the application. DETAILED DESCRIPTION

[0154] The application will be described in further detail below with specific embodiments. The embodiments of the application are given for illustrative and descriptive purposes only, and are not exhaustive or limiting of the application to the forms disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art. Embodiments are chosen and described in order to best explain the principles of the application and its practical application, and to enable others skilled in the art to understand the application for various embodiments with various modifications as are suited to the particular use contemplated.

[0155] Embodiment 1

[0156] Please refer to Figures 1-5 In the embodiment, a goaf grouting surface subsidence precision monitoring method is provided, characterized in that it comprises the following steps:

[0157] I. Establish a ground subsidence monitoring control network; research on grouting filling goaf surface precision leveling observation technology Establish a ground subsidence monitoring control system, lay out a high-precision monitoring control network according to the second-order leveling measurement accuracy, monitor the surface residual subsidence of the goaf under the grouting condition according to the third-order leveling measurement accuracy; adopt overall observation of measuring points and key point review, etc. to build a point-surface combined subsidence monitoring network, ensure the field observation accuracy, and improve the data acquisition efficiency; establish an air-space monitoring network system; it integrates satellite remote sensing, low-altitude unmanned aerial vehicle and laser Lidar technology to collect data.

[0158] II. Observation data processing and analysis; the ground subsidence condition and the surface land use and cover change condition of the research area are obtained by means of InSAR technology, laser radar point cloud, optical remote sensing image processing technology and other technical means, the adjustment model is optimized, the field data processing is carried out, the surface movement and deformation characteristics of the grouting influence area are analyzed, and the measured basis is provided for the research on the surface residual deformation law of the grouting filling goaf.

[0159] III. Similar material simulation experiment: the overburden strata movement law of the grouting filling goaf is researched by means of similar material simulation experiment, numerical simulation calculation and mechanical analysis, and the overburden strata movement mechanism of the grouting filling goaf is revealed

[0160] Four, experimental analysis: experimental test of soil, rock, filling body and filling body and surrounding rock composite physical and mechanical properties, obtain the mechanical properties of rock in goaf under different confining pressure and unloading rate and parameters, analyze and evaluate the mechanical properties of overburden before and after grouting and filling in goaf; Take the caving zone rock and the filling body in the goaf during construction, carry out hierarchical loading creep test, obtain the creep curve and long-term strength change rule, establish the corresponding constitutive model and carry out parameter identification;

[0161] Five, through the processing and analysis of observation data and simulation experiment data, the coupling relationship between residual deformation of goaf and compression deformation of building foundation is analyzed.

[0162] The high-precision ground subsidence monitoring system adopts the second-order leveling measurement precision to layout the high-precision control network, and monitors the ground residual subsidence in the key research area under the condition of goaf grouting according to the third-order leveling measurement precision.

[0163] (1) According to the actual situation of the survey area, a total of second-order closed leveling routes are arranged on the ground, and the distance between adjacent measuring points is 136m-635m.

[0164] (2) According to the actual situation of the survey area, a total of third-order leveling routes are arranged on the ground, and the distance between adjacent measuring points is 40m-60m.

[0165] In step one, the space-air monitoring network system is based on the space-air integrated monitoring technology system of Beidou, which integrates the GNSS and InSAR technology of precise and accurate ground deformation monitoring technology. The advantages of GNSS (Beidou high-precision positioning) and InSAR technology are integrated to realize the precise and accurate real-time monitoring of the ground deformation in the grouting engineering area, and to obtain the three-dimensional deformation field of the research area, which provides data support for analyzing the ground deformation law of the grouting engineering area under the influence of multiple factors; The regional fine ground model of unmanned aerial vehicle+LiDAR is used to construct continuous high-precision DEM and DSM models by using unmanned aerial vehicle carrying high-precision LiDAR, to directly describe and express the terrain surface, and to obtain the overall deformation result and deformation trend of the monitoring area; The ground cover change and historical evolution law of ground subsidence are monitored by using remote sensing technology, based on multi-temporal high-resolution remote sensing data, the historical evolution data of land, vegetation and other ground cover in the research area are accurately extracted, and the evolution law is obtained; The average deformation rate information and deformation history information of the ground target in the research area are obtained by using SBAS-InSAR and other data processing technologies.

[0166] The work adopts the method of combining D-InSAR with coherent target time series analysis to process the selected satellite radar data. The InSAR data processing process includes data collection, radar data registration, D-InSAR data processing and effective data selection, data quality evaluation, etc.

[0167] 1)Data collection

[0168] InSAR data, DEM data and basic geographic information data for mapping analysis were collected in the work.

[0169] (1) The radar data used in the InSAR ground subsidence information extraction work is sourced from the RadarSAT-2 satellite, and the precise orbit parameters of each period of data are obtained.

[0170] (2) The DEM data related to the work area in the InSAR data processing work was collected, and the DEM data obtained by SRTM was used, with a horizontal resolution of 30 meters, meeting the requirements of terrain compensation in InSAR processing.

[0171] (3) Other optical remote sensing image data, basic geographic information data, etc. for assisting ground subsidence status analysis, result mapping, etc.

[0172] 2) Radar data registration

[0173] Obtain the precise orbit data of each period of data, use orthogonal correlation technology to achieve precise registration between each primary and secondary image, and the registration accuracy of the whole scene data is required to be higher than 0.2 pixels. Register the DEM data with the RadarSAT-2 primary image as the reference, remove the terrain phase in subsequent interferometric measurement through DEM data information and radar imaging geometry information. Obtain the SRTM DEM data covering the test area, and according to the SAR image imaging geometry and other information, perform coordinate system conversion of external DEM, i.e. geocoding, to establish the conversion lookup table of the coordinate system between SAR image and external DEM. Through the texture matching of DEM simulated SAR image and SAR image, the matching accuracy is higher than 0.2 pixels, and the lookup table is corrected. Through the lookup table, the DEM can be sampled to the SAR radar imaging coordinate system, and then combined with the satellite orbit information to simulate the correct terrain phase, complete the difference processing of the interferometric phase. Similarly, in the later result output, the deformation information in the SAR imaging coordinate system can be sampled to the DEM geographic coordinate system, to complete the analysis and mapping work of the later results.

[0174] 3) Selection of interferometric image pairs and D-InSAR processing

[0175] Precise orbit data is added to each period of data to reduce the noise and terrain error caused by orbit error.

[0176] (1) Selection of interferometric image pairs

[0177] According to the baseline distribution characteristics of the acquired SAR data sequence, combined with the basic characteristics of regional ground deformation, each data and its adjacent 5 data form an interferometric image pair sequence, and 180 interferometric image pairs are formed from 39 data. According to the time sampling of 24-day data, the time baseline is between 24-120 days. Due to the missing of the middle part of the data, the time baseline of part of the interferometric pairs is longer, and is less than 240 days.

[0178] (2) Interference fringe sequence generation

[0179] The interferometric image pairs meeting the baseline requirements are subjected to interference processing to generate an interference fringe sequence. According to the characteristics of the noise phase, the interference fringe is subjected to adaptive filtering processing.

[0180] (3) Differential interference processing

[0181] Precise orbit parameters are introduced for baseline estimation, flat ground is removed according to the geometric imaging model, DEM sampled to radar coordinate system is used, and terrain phase is simulated through the geometric imaging model to obtain a differential interferogram containing deformation phase, atmospheric phase and noise phase.

[0182] (4) Interference fringe screening

[0183] The differential interference fringe is screened to exclude the interferograms with obvious strong convection influence and the interferograms with obvious residual interference phase stripes.

[0184] 4) Coherent target time series analysis

[0185] The interference phase of the coherent point target is subjected to time series analysis, the atmospheric fluctuation, DEM error and noise are estimated according to the space-time characteristics of each phase component, and are separated from the differential interference phase one by one, and finally the deformation rate of each coherent point target is obtained. The accuracy of the annual deformation rate obtained can reach millimeter level.

[0186] The differential interference phase of the coherent point target is extracted by using the differential interference phase obtained in the early stage, the spatial phase unwrapping is used for the point target phase, the Delaunay triangular network is constructed according to the positions of the points, the residual points in the triangular network are judged, the positive and negative residual point pairs are connected by using the minimum cost flow algorithm, and the unwrapping result of each point is obtained.

[0187] For the differential interference phase sequence of the coherent point target, two-dimensional periodic function is used for parameter inversion to estimate the adjacent point

[0188] The moving speed and height error correction are carried out through multiple error correction iterations, the correlation coefficient of the two-dimensional periodic function is 0.8-0.85, and the linear deformation rate is accurately extracted.

[0189] 5) Data comprehensive analysis and mapping

[0190] (1) Coordinate look-up table created by DEM registration, sampling the above point target deformation rate information to the common coordinate system.

[0191] (2) Comprehensive application of image processing and GIS technology to analyze InSAR ground deformation information and prepare InSAR monitoring ground subsidence map.

[0192] 6) Quality control of data processing

[0193] InSAR processing is very strong, and each step is checked by 100% human-computer interaction and mutual inspection between operators, mainly including:

[0194] (1) Raw data quality

[0195] Precise monitoring of surface movement and residual deformation law of goaf grouting and filling

[0196] Comprehensive inspection of radar data obtained from each scene, including data type, coverage, phase, beam mode, etc. The data is read into professional processing software, and the correct reading, no bad points in the map, and no loss of phase information.

[0197] (2) Image registration (radar data and DEM data)

[0198] Select the same main image to register radar data, and the registration accuracy distance and direction are higher than 0.2 pixels. The registration accuracy of DEM and main image reaches 0.2 pixels. If the registration accuracy does not meet the requirements, it cannot form effective phase in the subsequent differential interference processing.

[0199] (3) Differential interference processing

[0200] Differential interference processing includes interference of two SAR data, terrain phase removal, flat phase removal, and individual inspection of the formed differential interference phase diagram to ensure the accuracy of differential coherence.

[0201] (4) Residual error correction

[0202] Residual terrain phase removal, orbit error removal, and exclusion of significant atmospheric error interference phase for the differential phase diagram formed after differential interference processing.

[0203] (5) Selection of coherent target points

[0204] Check the accuracy of the position of the coherent point target and whether the coherent point target density meets the subsequent time analysis requirements. The point density in urban areas is required to be more than 100 per square kilometer.

[0205] (6) Coherent target time series analysis

[0206] Carry out coherent point target time series analysis to extract land subsidence rate information. Check the phase unwrapping generated by time series analysis

[0207] File, residual file and deformation file, ensure that each layer is unwrapped correctly, there is no discontinuous phase, and the residual file has no blocking phenomenon.

[0208] (7) Geocoding

[0209] Geocode the deformation information and check the correctness of the file's geographic location.

[0210] (8) Format conversion

[0211] Convert the deformation information to shp data format for analysis and mapping, and check the correctness of shp geometry and attributes.

[0212] Step one, field data collection

[0213] Use the vertical take-off and landing fixed-wing unmanned aerial vehicle CW-25E to carry LiDAR system for aerial photogrammetry, and the ground station is GCS-202

[0214] Step two, point cloud data preprocessing

[0215] In data processing, the profile method is used to check data quality, and the specific steps are as follows:

[0216] 1) Render the point cloud into different colors according to the elevation value, that is, point cloud coloring. Point cloud coloring is to give the original LiDAR point cloud data without RGB color information with actual texture color. The ground information in the survey area can be displayed intuitively in the point cloud view, and the RGB color information can be used to assist point cloud classification in the later period. There are mainly two forms to realize point cloud coloring;

[0217] (1) Color according to the original image taken in the survey area, camera calibration file, POS file after aerial triangulation and color range;

[0218] (2) Color according to the orthophoto generated after aerial triangulation or fast mosaic. Data processing uses the method of implementing point cloud coloring based on orthophoto.

[0219] 2) Use the profile function to make a profile in the main interface of each overlapping area of the flight strip, and observe the profile view to see if there is obvious layering. If it is found that the point cloud is layered, it needs to be reprocessed for point cloud calculation and flight strip adjustment. In addition, the quality inspection function in the intelligent laser module of the unmanned aerial vehicle manager can also be used to check the height difference of the point cloud between the flight strips, which can also reflect whether the point cloud is layered from the side.

[0220] 3) Coordinate transformation

[0221] The required point cloud coordinate system is CGCS2000 national geodetic coordinate system, Gauss 3-degree zone projection, central meridian 117°, 1985 national height datum. After point cloud adjustment and color assignment, coordinate conversion is carried out, mainly involving parameter calculation and coordinate conversion. The residuals dN, dE and dU of the seven parameters to be solved must be controlled within the limit error, and the detailed steps are as follows:

[0222] (1) The initial point cloud is in the WGS84 system, and now the point cloud in the CGCS2000 national geodetic coordinate system is required. Parameter conversion between the two coordinate systems is required. Three parameter points are measured in the field, and three pairs of coordinates of the three points in the WGS84 coordinate system and the CGCS2000 national geodetic coordinate system are known. The seven parameters of the two coordinate systems are calculated, and the residual information dN, dE and dU is checked to see if it meets the requirements.

[0223] (2) The coordinate conversion parameter configuration of the survey area is carried out, the Bursa seven parameters obtained by parameter solving are added, the conversion parameter configuration is completed, and the point cloud is converted from the source coordinate system (WGS84 coordinate system) to the target coordinate system (CGCS2000 national geodetic coordinate system).

[0224] 4) Standard point cloud data output

[0225] After the above preprocessing operation of unmanned aerial LiDAR point cloud data, point cloud data for production can be generated, and the obtained point cloud results are generally stored in the form of “*.las” format. Using ground point cloud data and image data, three-dimensional models, DEM, DSM and DOM digital products can be generated as needed.

[0226] Step two, point cloud data processing

[0227] After the above point cloud data preprocessing, standard point cloud data in “.las” format can be obtained. Unmanned aerial LiDAR is affected by topographic environmental features, equipment accuracy, characteristics of ground objects and other factors during operation. The obtained point cloud data contains certain noise points, which will affect the accuracy of subsequent data processing and monitoring results. In addition, unmanned aerial LiDAR scanning is scanning all ground objects such as buildings, vegetation, power facilities, etc. Point cloud data contains surface information of all ground objects. In order to fully utilize point cloud data and achieve better

[0228] The application effect of goaf grouting filling ground surface movement precision monitoring and residual deformation law research needs to classify the obtained point cloud data according to actual needs.

[0229] 1) Point cloud denoising

[0230] The acquired mine standard format las point cloud can be visualized, and a large number of noise points can be obtained. Noise points exist above and below the ground point cloud. In order to ensure the accuracy of the DEM data product generated by using the point cloud data, it is necessary to remove the noise point group.

[0231] The point cloud denoising function in the point cloud module of Feima unmanned aerial vehicle manager is selected to denoise the point cloud data of the entire project area.

[0232] 2) Point cloud filtering and classification

[0233] Point cloud filtering refers to the process of removing non-ground points (vegetation, buildings, bridges, power facilities, etc.) in point cloud data in order to extract a digital elevation model from point cloud data. This process is called filtering of laser radar data.

[0234] Point cloud classification refers to the process of categorizing non-ground points, which is used for feature extraction and three-dimensional reconstruction of buildings.

[0235] The purpose of the research on airborne laser radar data processing is to monitor the ground settlement before and after grouting in the project area. Therefore, the main goal of laser point cloud filtering and classification is to distinguish between ground points and non-ground points. The most basic principle of filtering is that there is a height mutation between the laser foot point and the surrounding point cloud. For example, there is a significant difference in height between the laser points corresponding to buildings, streetlights, utility poles, vegetation, and the laser points corresponding to adjacent features. The height of different parts of a tree is also different. Point cloud classification includes automatic classification based on classification algorithms and manual classification based on manual editing.

[0236] Automatic classification of point cloud is the use of algorithms or algorithm combinations based on the reflectivity, echo frequency, shape characteristics, etc. of different ground objects to automatically classify point clouds representing different types of features, mainly separating into ground points, non-ground points, and air points, low points, and noise points. After automatic classification, the point cloud may have misclassification, under-classification, or over-classification. Manual editing is required to reclassify points in the filtered point cloud that do not belong to the ground, and the classified points are stored in the corresponding category through classification editing. The Feima unmanned aerial vehicle manager intelligent point cloud module is based on the irregular triangular grid encryption method for ground point classification. The irregular triangular grid encryption method is to divide the point cloud of the entire study area according to the terrain slope, select the lowest point in each block to generate a sparse triangular grid, and then encrypt layer by layer until all ground points in the block are classified. Finally, all classified single-block point clouds are merged. The specific operation process is to set the point cloud data to editable mode, select the ground point option in the ground point tab in the data editing menu bar, and perform automatic classification of the ground point cloud. During classification, parameters such as maximum building size, maximum slope threshold, maximum iteration angle, and maximum iteration distance need to be set according to the terrain characteristics of the project area.

[0237] The ground point cloud extracted after automatic classification of point cloud, and the comparison of the results before and after automatic classification of building area ground point cloud (partial area). It can be seen that the automatic classification of ground point cloud has obvious effect, and the point cloud has strong penetration ability through vegetation. Many points penetrate through the vegetation to reach the ground. Through point cloud automatic classification, the vegetation and buildings on the ground can be well removed from the point cloud

[0238] 3) Digital product generation

[0239] After point cloud denoising and automatic classification of standard point cloud by using the intelligent point cloud module of unmanned aerial vehicle manager, high-precision laser point cloud data for producing digital products is obtained. By using the data editing function in the intelligent point cloud module, the point cloud after denoising processing can be generated DSM, and the ground point cloud separated by automatic classification can be further used to generate DEM. For the area with poor automatic classification effect, manual classification can be performed through regional smoothing, online classification and offline classification. DSM and DEM digital products are generated by using the point cloud obtained after denoising and classification processing. After point cloud denoising and automatic classification, DSM is generated by using the obtained ground point cloud. There is no non-ground point cloud such as ground buildings and vegetation in DSM, and the classification effect is significant.

[0240] However, the ground point cloud extracted by automatic classification contains some non-ground point cloud. Some point clouds are obviously lower or higher than the ground. The non-ground point cloud that cannot be removed by point cloud automatic classification can be manually removed by artificial classification method. In the intelligent point cloud module of Feima unmanned aerial vehicle manager, various manual classification methods such as brush classification, offline classification, online classification and smoothing classification are given. In profile mode, a polyline is manually drawn as a reference. The software will automatically classify the point cloud below or above the polyline as target point cloud. The data generated after the data processing of the target point cloud mainly includes LiDAR point cloud data, DOM, DEM and DSM.

[0241] 4) Data accuracy check mainly includes: point cloud data density and point cloud data plane and elevation accuracy check.

[0242] By using the point density measurement function in the intelligent laser module of Feima unmanned aerial vehicle manager, the point density of point cloud data can be measured. The measurement box side length value can be set to measure the point density in different areas. Three types of areas in the project area are selected, and the point density is measured. The three types of areas are building dense area, vegetation dense area and ground area. The accuracy check function can check the elevation accuracy of point cloud. By importing the check point file and comparing it with the point cloud data obtained by point cloud calculation, the elevation difference between the point cloud data and the check point coordinates is obtained, and the accuracy check of point cloud data is realized.

[0243] In step five, the coupling relationship between the residual deformation of the goaf and the compression deformation of the building foundation is analyzed, and the detailed steps are as follows:

[0244] 1) Residual surface deformation calculation method under the condition of goaf grouting

[0245] After coal mining, any point in the subsidence basin will undergo a long duration of movement process, and the movement period can be roughly divided into initial period, active period and recession period. After the end of the recession period, the surface enters the subsidence stable period, at this time the surface is still not completely stable, because there are a large number of cavities, parting and cracks in the stratum, even if the grouting and filling is carried out, it is still impossible to completely avoid the existence of underdense space, and the filling material will also produce dry shrinkage deformation after dehydration, under the action of overburden load, the surface will still produce a certain residual subsidence deformation.

[0246] Research shows that the form of residual deformation of goaf is mainly the compression deformation of overlying strata. The rock mass in caving zone is in a state of loose and fragmented structure, with large interstitial space and good connectivity, and the density of broken rock mass in goaf can be improved by grouting and filling,

[0247] The filled body after solidification forms a cemented composite with broken rock and surrounding rock, and under the influence of overburden load and dehydration, the deformation characteristics of the cemented composite are compression and displacement under external force. The rock mass in the fracture zone is in a state of block and layer structure, with large block size, high stiffness and strength; above the boundary of the goaf, the rock blocks interlock to form a masonry beam type semi-arch structure, and due to the difference in fracture block size, layer thickness and stiffness, there are a large number of parting and cracks between layers, which can be filled by grouting and filling to densify the parting and cracks in the rock layer, improve the overall performance and bearing capacity of the rock layer, and slow down the movement and deformation of the rock layer, and the deformation characteristics are mainly the compression of rock block structure and cracks. The rock mass in the bending zone is in a relatively complete layer structure, and the lower part often has parting phenomenon, which can be further controlled after grouting and filling of the parting, and the deformation characteristics are the compression of parting and cracks under the action of overburden pressure. Grouting and filling in goaf can slow down the residual deformation of goaf surface, and the form of surface movement is relatively flat residual subsidence, which is in the shape of a basin with large middle and small edge above the goaf, and the subsidence basin still conforms to the normal curve distribution rule, and the probability integral method mathematical model can be used for residual surface deformation calculation.

[0248] Principle of probability integral method, the probability integral method is an analytical method of random medium theory, which originated in 1954 from the mining subsidence random medium theory proposed by Polish scholar Litvinisn. This method generalizes the rock mass as a random medium composed of a large number of loose granular substances, and simplifies the movement of rock strata and surface in goaf as the random movement of a large number of granular media, from the statistical point of view, the entire mining under any mining condition is decomposed into a number of micro-unit mining, and the influence of each micro-unit mining on the movement of rock strata and surface in goaf is equivalent to the influence of the whole mining on the rock strata and surface.

[0249] (1)Establish the probability integral method theoretical model, in the theoretical model, the medium particles are assumed to be some small balls with the same size and uniform mass, and are packed in the same size and uniform arrangement of the grid, as shown in Figure 6 If the small ball in the lower grid is removed, the probability of the small ball in the upper two adjacent grids rolling into this grid is 1 / 2 due to the action of gravity, and thus the probability distribution diagram of the particle movement is obtained by analogy, as shown in Figure 7 If the grid and the particle have no goaf grouting filling ground surface movement precision monitoring and residual deformation law research limit small, the medium movement probability distribution tends to a smooth curve, the probability density curve, as shown in Figure 8

[0250] Through mathematical derivation, the expression of the subsidence value We(x,z) of the point (x,z) caused by unit mining can be obtained.

[0251]

[0252] In the formula:

[0253] W e (x,z)——the subsidence value caused by unit mining, unit mm;

[0254] r z ——the main influence radius, unit m,

[0255]

[0256] For the ground surface, z is equal to the mining depth, and A is a constant representing the size of the medium unit, so that r z is a constant, and let r z be a constant r, then the expression of the subsidence basin of the ground surface unit is:

[0257]

[0258] (2)Residual ground surface deformation prediction of goaf grouting and filling

[0259] Apply the probability integral method formula, take the mining thickness of a point (s,t) in the calculation area of the goaf as Md(s,t), and the spatial position relationship between the calculation area and the ground surface is shown in Figure 2 , then the subsidence of any point (x,y) on the ground surface is:

[0260]

[0261] In the formula:

[0262] l A , l B ​Length and width of the goaf calculation area, unit: m;

[0263] M d Mining thickness of coal seam, unit: m;

[0264] q c Residual surface subsidence coefficient of goaf grouting and filling;

[0265] α - Coal seam dip angle;

[0266] r - Main influence radius, unit: m.

[0267] Then the tilt, curvature, horizontal movement and horizontal deformation of a given φ direction on the surface are:

[0268]

[0269]

[0270] (3) Analysis of residual surface subsidence calculation parameters

[0271] The residual surface movement and deformation of goaf grouting and filling can be predicted by probability integral method. Generally speaking, the prediction results of corresponding parameters should be reliable by fully verified by measured data. However, due to the long duration and small value of residual surface subsidence of goaf grouting and filling, it is generally difficult to grasp its overall development law by measured method. Even if there are observation points on the surface of goaf, after the mining movement is basically stable, the observation is continued for several years,

[0272] When the settlement of the observation point is monitored to be within the error range of the observation instrument, in fact, at this time, it cannot be determined that the residual surface subsidence of the goaf has ended. When the goaf is "activated" under the action of water filling or water loss change, near mining influence, earthquake or ground stress, the surface of the goaf will again subsidence, and the time of these conditions or actions is difficult to predict, so the residual surface deformation has unpredictability. Considering the particularity and complexity of the surface movement law of goaf grouting and filling, as well as the subsidence reduction mechanism of goaf grouting and filling, the residual surface deformation of goaf grouting and filling can be calculated by the way of calculating the residual subsidence coefficient according to the filling rate of goaf.

[0273] (a) Residual surface subsidence coefficient qc of goaf grouting and filling

[0274] The residual subsidence coefficient is mainly related to the subsidence coefficient in the surface movement duration and the duration after the end of mining. The greater the subsidence coefficient in the surface movement duration, the better the compaction of the overlying strata during the mining process, and the smaller the residual surface deformation after the surface is stable, and vice versa. The longer the working face stop mining time, the smaller the residual subsidence.

[0275]

[0276] q - the subsidence coefficient in the surface movement duration;

[0277] q1 - the residual surface subsidence coefficient under general mining conditions;

[0278] k - adjustment coefficient, generally 0.5-1.0;

[0279] t - time from the end of mining;

[0280] The residual deformation of the longwall fully-mined working face can be regarded as the compression deformation of the broken roof rock mass of the goaf under the pressure of the overlying strata. The size of the compression deformation is related to the compression rate of the broken roof rock mass. The purpose of goaf grouting is to fill the voids of the goaf and the rock voids in the caving zone, increase the compaction of the rock in the caving zone, and at the same time fill the joints, fissures, and bedding in the fractured zone to improve the overall strength of the rock strata. Therefore, the surface movement deformation under the condition of goaf grouting and filling can be calculated by the goaf compaction rate

[0281] The residual subsidence coefficient can be calculated by the inverse calculation method, and the residual surface subsidence coefficient under the condition of goaf grouting and filling can be calculated by the following formula:

[0282]

[0283] In the formula: q c - residual surface subsidence coefficient under the condition of goaf grouting and filling;

[0284] ρ c - goaf filling rate;

[0285] ρ y - compaction rate of the goaf filling body.

[0286] The compaction rate of the goaf filling body is related to the mining depth of the coal seam, the lithology of the overlying strata, the filling material, the slurry ratio, and the construction technology. The greater the mining depth, the greater the pressure acting on the filling body, and the greater the compression strain of the filling body;

[0287] The softer the overburden is, the greater the pressure on the filling body is, and the greater the compression of the filling body is. The compaction rate can be measured by compression test, and the calculation method can be calculated according to formula 11, generally taking about 0.1

[0288]

[0289] In the formula:

[0290] V - the volume of the filling body filled in the goaf, unit m 3 ;

[0291] V' - the volume of the filling body after compression, unit m 3

[0292] (b) inflection point offset coefficient S / H

[0293] After the end of mining in the goaf, there are a large number of cavities, cracks and roof cantilever structures near the goaf boundary, which causes the inflection point of the subsidence curve at the time of the surface movement to be stable to be located inside the goaf at a certain distance, which is the inflection point offset, represented by the inflection point offset coefficient S / H. According to the "three-under" coal mining specification, when the overburden type is hard, medium-hard and soft, respectively, S / H at the time of the surface basic stability is 0.31-0.43, 0.08-0.30, and 0-0.07, respectively.

[0294] Precise monitoring of surface movement and residual deformation law of goaf grouting and filling

[0295] Through field engineering practice, it is found that the goaf under the project area has been compacted for many years, and the cavities and cracks near the goaf boundary are more developed. Through goaf grouting and filling, the cavities and cracks near the goaf boundary can be effectively filled, which can support the cantilever structure of the roof above to a certain extent and prevent the roof from continuing to collapse. Therefore, the inflection point offset coefficient when calculating the residual surface subsidence of the goaf can be analogized to the parameters at the time of the surface basic stability.

[0296] (c) main influence angle tangent tanβ

[0297] The size of the main influence angle tangent tanβ is related to the mining method, roof management method, overburden lithology and thickness, etc. After the mining method and roof management method are determined, tanβ is mainly related to the lithology and thickness of the overburden. Goaf grouting and filling will not significantly change the physical and mechanical properties of the overburden, therefore, the main influence angle tangent tanβ when calculating the residual subsidence will not change significantly with the extension of time after the end of mining, and its value can be referred to the tanβ of conventional subsidence deformation. According to the "three-under" coal mining specification, when the overburden type is hard, medium-hard and soft, respectively, tanβ at the time of the basic stability of mining is 1.2-1.91, 1.92-2.40, and 2.41-3.54, respectively.

[0298] (d) Horizontal movement coefficient b

[0299] The horizontal movement coefficient b is mainly related to the properties of overburden, and has little to do with other factors. The horizontal movement coefficient b in the calculation of residual surface subsidence can refer to the value of conventional subsidence deformation. According to the "three-under" coal mining specification, the value range of horizontal movement coefficient is 0.2-0.4, generally 0.3.

[0300] (e) Mining influence propagation angle θ

[0301] The mining influence propagation angle θ depends on the coal seam dip angle and the mining influence propagation coefficient K, which is mainly related to the overburden lithology, coal seam mining depth, etc. The influence of gob grouting and filling on it is limited. The mining influence propagation angle θ in the calculation of residual surface subsidence can refer to the value of conventional subsidence deformation θ. According to the "three-under" coal mining specification, θ = 90°-kα (α is the coal seam dip angle), and the value of k is 0.7-0.8, 0.6-0.7, and 0.5-0.6 when the overburden type is hard, medium hard, and weak, respectively.

[0302] Building load under the compression settlement calculation model of foundation soil

[0303] When engineering construction is carried out above the goaf, the newly built building (structure) will be affected by the residual deformation of the goaf and the compression settlement deformation of the foundation soil under the action of building load, therefore, the residual surface deformation above the goaf should include the superposition of the residual movement deformation of the goaf and the compression settlement of the soil under the action of load.

[0304] (1) Spatial stress analysis of foundation soil under load

[0305] The influence range of building rectangular load (length-width a×b) on the underlying foundation soil is approximately an ellipsoid, for the convenience of calculation, the maximum influence range of soil can be regarded as a cuboid with length-width-height la, lb, Hz, see Figure 3where la and lb are the length and width of the soil body affected by the load, and Hz is the load-affected depth. According to the Regulations for Precision Monitoring of Ground Surface Movement and Residual Deformation Law of Grouting Filling in Goaf, when the additional stress caused by the building load in the foundation is equal to 20% of the self-weight stress of the foundation rock-soil layer at the corresponding depth, the influence of the additional stress on the soil at the depth can be ignored, but when there is high compressibility soil or a goaf collapse, fracture zone or other instability factors below, the calculation should be performed to the position where the additional stress is 10% of the self-weight stress of the foundation rock-soil layer, and the depth is the building load-affected depth minus the foundation burial depth. Since the additional stress is the largest at the center point of the rectangular foundation, the load-affected depth is the largest below the center point, and the additional stress gradually decreases from the center of the rectangle to the outside of the foundation, until the influence of the additional stress on the soil can be ignored in a certain range of soil outside the foundation, and the influence depth is 0,

[0306] Therefore, the critical point connecting the load-affected depth of 0 calculated on both sides of the ground surface along the center lines of the long and short sides of the foundation is taken as the length of la and lb, and Hz takes the maximum value of the load-affected depth, i.e. the load-affected depth below the center point of the rectangle. The spatial relationship among la, lb and Hz is shown in Figure 4 . In the affected space range of the soil, take a point M with coordinates (x, y, z), and the rectangular foundation is subjected to uniform load P. At point M, a unit body is taken, and the additional stresses in the vertical and horizontal directions at the unit body can be calculated according to formula 12

[0307]

[0308] where σ fz(x,y,z) , σ fx(x,y,z) , σ fy(x,y,z) are the additional stresses in the x, y and z directions at point (x, y, z).

[0309] σ z(x,y,z) , σ x(x,y,z) , σ y(x,y,z) are the additional stress coefficients in the x, y and z directions at point (x, y, z).

[0310] (2) Layered settlement model of foundation soil

[0311] The foundation soil has the characteristics of natural soil layering, and the compression characteristics of the soil are considered. The layered summation method in soil mechanics is applied to calculate the foundation settlement. The model has the following assumptions:

[0312] (a) Each layer in the foundation is vertically compressed without lateral expansion.

[0313] (b) The thickness of each layer should not be too large, generally 1-2 m, and the stress from the top to the bottom of each layer is changing, and the average value of the stress at the top and bottom is approximately taken in the calculation.

[0314] (3) Calculate the soil deformation in the load influence range.

[0315] The settlement of a point (x, y, 0) on the surface of the foundation soil is the cumulative compression settlement of each unit soil from the point to the depth of z = Hz. The settlement value of the unit at point M (x, y, z) can be calculated by the following formula 13 for dS (X,Y,Z) Integrating on z = [0, Hz], the settlement at point (x, y, 0) can be calculated by the following formula 14:

[0316]

[0317] In the formula:

[0318] dS (X,Y,Z) — The vertical settlement of the unit soil at point (x, y, z);

[0319] E S(X,Y,Z) — The compression modulus of the unit soil at point (x, y, z);

[0320] S (x,y,0) — The settlement at point (x, y) on the surface (z = 0) of the foundation soil.

[0321] Since the foundation soil has the characteristics of natural soil layering, from the z = 0 plane to the z = Hz plane, it can be divided into m1 soil layers, then S (x,y,0) It can also be calculated by formula 14:

[0322]

[0323] In the formula:

[0324] — The average additional stress of the i-th soil layer below point (x, y, 0);

[0325] — The average additional stress coefficient of the i-th soil layer below point (x, y, 0);

[0326] Δz — The thickness of the i-th soil layer below point (x, y, 0);

[0327] E s — The compression modulus of the i-th soil layer below point (x, y, 0);

[0328] m1 — The number of soil layers in the load influence depth range;

[0329] p — The average pressure of the foundation, unit kg / m2.

[0330] The settlement of the building foundation under the action of load can be considered as the superposition of the compression settlement of the foundation soil and the residual surface subsidence, i.e. S总 = S d + S w , the origin of the selected coordinate system is different, see Figure 5 , the origin of the selected coordinate system is different, see

[0331]

[0332] The compression of the soil body under the action of the building load is mainly vertical settlement, and the influence on lateral deformation is very small and can be ignored. The deformation of the foundation soil layer under the load of the upper building causes the settlement of the building foundation. When the site soil is solid, the settlement of the foundation is small, and the influence on the construction project is slight; but if the foundation is soft soil layer and the thickness is uneven, or the load of the upper structure varies greatly, the foundation will have obvious uneven settlement, which will cause adverse effects on the normal use of the building. The settlement and horizontal movement of the foundation soil body are mainly the cooperative settlement of the residual deformation of the goaf and the compression deformation of the foundation soil body, so the inclination, horizontal movement, curvature and horizontal deformation of a given p direction on the surface of the foundation soil body can be represented by the following formula

[0333]

[0334] (4) Building foundation settlement deformation calculation under the condition of goaf grouting filling

[0335] The settlement of the building foundation under the condition of goaf grouting filling mainly includes two parts of residual surface subsidence of goaf and compression settlement of foundation soil body, and the specific calculation is shown in formulas 17-20. In the calculation, the residual surface subsidence coefficient can be calculated according to formula 9, or it can be selected according to experience, and other parameters can be similar to the residual surface movement parameters of longwall mining in the "three-under" coal mining specification

[0336] The compression settlement of the foundation soil body is calculated by the layer summation method. The ground building (structure) is mostly a cuboid, and its foundation soil body stress has symmetry, so the settlement of 1 / 4 area of the foundation soil body cuboid space can be calculated. The settlement of a point on the surface of the foundation soil body is the cumulative compression settlement of each unit soil body on the straight line from the point to the load influence depth (the additional stress is 5% of the weight of the soil layer as the criterion for determining the influence depth). The soil body needs to be layered according to the type and thickness of the soil body, and the layering thickness should not be too large (not more than 3m), and the settlement of an arbitrary point on the surface of the soil body is calculated by formula 14.

[0337] Obviously, the described embodiments are only some of the embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art and related fields without creative labor should belong to the protection scope of the present application. The structures, devices and operation methods not specifically described and explained in the present application are implemented according to the conventional means in the art, if not specifically described and limited.

Claims

1. A method for precise monitoring of surface settlement during grouting in goaf areas, characterized in that, Includes the following steps: I. Establish a ground subsidence monitoring and control network; research on precise leveling observation technology for grouting and filling goaf areas; establish a ground subsidence monitoring and control system, deploy a high-precision monitoring and control network with second-order leveling accuracy, and monitor residual surface subsidence in the goaf area under grouting conditions with third-order leveling accuracy; adopt methods such as comprehensive observation of measuring points combined with verification of key measuring points to construct a subsidence monitoring network combining points and areas, ensuring the accuracy of field observations and improving data acquisition efficiency; establish an air-space monitoring network system; which integrates satellite remote sensing, low-altitude UAVs and laser LiDAR technology for data acquisition. II. Data Processing and Analysis: The study area's ground subsidence and land use cover changes were obtained through various technologies such as InSAR, lidar point cloud, and optical remote sensing image processing. The adjustment model was optimized, field data was processed, and the characteristics of surface movement and deformation in the grouting-affected area were analyzed, providing experimental evidence for the study of residual deformation patterns in grouting-filled goaf areas. III. Similar Material Simulation Experiments: Using similar material simulation experiments, numerical simulation calculations, and mechanical analysis, the movement law of the overlying strata in grouting-filled goaf areas will be studied to reveal the movement mechanism of the overlying strata in grouting-filled goaf areas. IV. Experimental Analysis: The physical and mechanical properties of soil, rock, backfill, and the composite of backfill and surrounding rock were tested to obtain the rock mechanical properties and parameters of the goaf under different confining pressures and unloading rates. The mechanical properties of the overlying rock before and after grouting and backfilling of the goaf were analyzed and evaluated. Graded loading creep tests were conducted on rock blocks from the caving zone and the backfill of the goaf during construction to obtain creep curves and long-term strength variation laws. Corresponding constitutive models were established and parameters were identified. V. By processing and analyzing observation data and simulation experimental data, the coupling relationship between residual deformation in the goaf and compressive deformation of the building foundation is analyzed.

2. The method for precise monitoring of surface settlement during grouting in goaf areas according to claim 1, characterized in that, In step one, the air-space monitoring network system is an integrated air-space monitoring technology system based on BeiDou. It integrates GNSS and InSAR technologies to conduct precise and accurate monitoring of surface deformation. This research aims to achieve precise and real-time monitoring of surface deformation in the grouting project area, acquiring the three-dimensional deformation field of the surface in the study area, and providing data support for analyzing the surface deformation patterns under the influence of multiple factors. A regional refined surface model using UAV + LiDAR is employed. Using a UAV equipped with high-precision LiDAR, continuous high-precision DEM and DSM models are constructed to directly describe and represent the terrain surface, acquiring the overall deformation results and deformation trends of the monitored area. Remote sensing technology is used to monitor the historical evolution of land cover changes and land subsidence. Based on multi-temporal high-resolution remote sensing data, historical evolution data of land cover such as land and vegetation in the study area are accurately extracted to obtain evolution patterns. Data processing technologies such as SBAS-InSAR are used to acquire the average deformation rate information and deformation history information of ground targets in the study area.

3. The method for precise monitoring of surface settlement during grouting in goaf areas according to claim 1, characterized in that, Step two employs a method combining D-InSAR and coherent target time series analysis to process the selected satellite radar data using InSAR. The InSAR data processing includes data collection, radar data registration, D-InSAR data processing and selection of valid data pairs, and data quality evaluation. Specifically, it includes the following steps: 1) Data collection: Essential InSAR data, DEM data, and basic geographic information data for cartographic analysis were collected during the work. 2) Radar Data Registration: Precise orbital data from each period is acquired. Orthogonal correlation technology is used to achieve precise registration between the main and secondary images, with a registration accuracy of over 0.2 pixels for the entire scene. DEM data is registered using the RadarSAT-2 main image as the reference. Topographic phase in subsequent interferometry is removed using DEM data information and radar imaging geometry. SRTM DEM data covering the test area is acquired. Based on SAR image imaging geometry and other information, coordinate system transformation (geocoding) is performed on the external DEM, establishing a coordinate system transformation lookup table between the SAR image and the external DEM. Texture matching between SAR images is simulated using the DEM, with a matching accuracy of over 0.2 pixels, and the lookup table is corrected. Using the lookup table, the DEM can be sampled into the SAR radar imaging coordinate system, and then combined with satellite orbital information to simulate the correct topographic phase, completing the differential processing of the interferometric phase. Similarly, in later output results, deformation information from the SAR imaging coordinate system can be sampled into the DEM geographic coordinate system to complete the analysis and mapping of the later results. 3) Selection of interferometric image pairs and D-InSAR processing: Add precise orbit data to the data of each period to reduce the impact of noise and terrain errors caused by orbit errors; 4) Coherent target time series analysis: Time series analysis is performed on the interferometric phase of coherent point targets. Based on the spatiotemporal characteristics of each phase component, atmospheric fluctuations, DEM errors and noise are estimated and separated from the differential interferometric phase one by one. Finally, the deformation rate of each coherent point target is obtained. The accuracy of the annual deformation rate can reach the millimeter level. Using the differential interferometric phase obtained in the early stage, the interferometric phase of the coherent point target is extracted. Spatial phase unwrapping is applied to the point target phase. A Delaunay triangulation is constructed based on the position of these points. The residual points in the triangulation are identified. The minimum cost flow algorithm is used to connect the positive and negative residual point pairs to obtain the unwrapping result of each point. For the differential interferometric phase sequence of coherent point targets, parameter inversion is performed using a two-dimensional periodic diagram to estimate the moving rate of adjacent points and the elevation error correction. Through multiple error correction iterations, the correlation coefficient of the two-dimensional periodic function is made to be 0.8-0.85, and the linear deformation rate is accurately extracted. 5) Data integration, analysis, and charting; 6) Quality control of data processing: InSAR processing is highly process-oriented, and each step is subject to 100% human-computer interaction inspection and mutual inspection among operators.

4. The method for precise monitoring of surface settlement during grouting in goaf areas according to claim 3, characterized in that, In step 1), data collection involves the following steps: (1) The radar data used for InSAR ground subsidence information extraction was obtained from RadarSAT-2 satellite, and the precise orbital parameters of each period of data were obtained. (2) The DEM data of the working area involved in the InSAR data processing was collected. The DEM data was obtained by SRTM with a horizontal resolution of 30 meters, which meets the terrain compensation requirements in InSAR processing. (3) Other optical remote sensing image data and basic geographic information data that assist in the analysis of ground subsidence and the compilation of results.

5. The method for precise monitoring of surface settlement during grouting in goaf areas according to claim 3, characterized in that, In step 3), the selection of interferometric image pairs and D-InSAR processing are detailed below: (1) Selection of interferometric image pairs: Based on the baseline distribution characteristics of the acquired SAR data sequence and combined with the basic characteristics of regional surface deformation, each data pair is composed of five adjacent data pairs. A total of 180 interferometric image pairs are formed from 39 data periods. Based on the sampling time of 24 days of data, the time baseline is between 24 and 120 days. Due to the missing data in the middle, the time baseline of some interferometric pairs is relatively long, all less than 240 days. (2) Interference pattern sequence generation: Interference processing is performed on the interferometric image pairs that meet the baseline requirements to generate an interference pattern sequence. Based on the characteristics of the noise phase, the interference pattern is subjected to adaptive filtering. (3) Differential Interferometry: Precise orbital parameters are introduced for baseline estimation. Flat areas are removed based on the geometric imaging model. The DEM sampled to the radar coordinate system is used to simulate the terrain phase through the geometric imaging model to obtain a differential interferogram containing deformation phase, atmospheric phase and noise phase. (4) Screening of interference patterns: Screen the differential interference patterns to exclude interference patterns that are significantly affected by strong convection and interference patterns with obvious residual interference phase fringes.

6. The method for precise monitoring of surface settlement in goaf grouting according to claim 3, characterized in that, In step 5), the specific steps for data synthesis, analysis, and charting are as follows: (1) Using the coordinate lookup table created during DEM registration, the deformation rate information of the above point targets is sampled into a commonly used coordinate system; (2) The InSAR surface deformation information was analyzed by comprehensively applying image processing and GIS technology, and InSAR monitoring ground subsidence maps were compiled.

7. The method for precise monitoring of surface settlement during grouting in goaf areas according to claim 3, characterized in that, The quality control of data processing in step 6) specifically includes the following: (1) Raw data quality: The study on the precise monitoring of surface movement and residual deformation of grouting and filling in the goaf area conducts a comprehensive check on the radar data of each scene, including data type, coverage, time phase, and beam mode. The data is read into professional processing software, and the data is read correctly, with no bad points in the map and no loss of phase information. (2) Image registration (radar data and DEM data): Select the same main image to register the radar data. The registration accuracy is higher than 0.2 pixels in both range and direction. The accuracy of DEM after registration with the main image reaches 0.2 pixels. If the registration accuracy requirement is not met, an effective phase cannot be formed in the subsequent differential interferometry processing. (3) Differential interferometry: Differential interferometry includes interferometry of two SAR data, terrain phase removal, and flatland phase removal. The resulting differential interferometric phase maps are checked one by one to ensure the accuracy of differential coherence. (4) Residual correction: The differential phase map formed after differential interferometry is subjected to residual topographic phase removal and orbital error removal, and the interferometric pairs with particularly significant atmospheric errors are excluded and no further calculation is performed; (5) Selection of coherent target points: Check the accuracy of the location of coherent target points and whether the density of coherent target points meets the requirements of subsequent time analysis. For urban areas, the point density requirement is more than 100 points / square kilometer. (6) Coherent target time series analysis: Conduct coherent point target time series analysis to extract ground subsidence rate information. Check the phase unwrapping file, residual file and deformation file generated by the time series analysis to ensure that each layer is correctly unwrapped, there are no discontinuous phases, and the residual file has no block phenomenon; (7) Geocoding: Geocoding the deformation information to check the correctness of the file's geographical location; (8) Format conversion: Convert deformation information into shapefile data format for analysis and plotting, and check the correctness of shapefile geometry and attributes.

8. The method for precise monitoring of surface settlement during grouting in goaf areas according to claim 1, characterized in that, In step two, the lidar point cloud data is preprocessed. The data processing employs a cross-sectional approach to check data quality. The specific steps are as follows: 1) Point clouds are rendered in different colors based on elevation values, i.e., point cloud colorization. Point cloud colorization assigns actual texture colors to raw LiDAR point cloud data that lacks RGB color information. The ground feature information in the survey area can be displayed intuitively in the point cloud view. Later, RGB color information can be used to assist in point cloud classification. There are two main forms of point cloud colorization: (1) Color assignment is performed based on the original images taken in the survey area, camera calibration files, POS files after aerial triangulation, and the color assignment range; (2) Coloring is performed on orthophotos or quick mosaics generated after aerial triangulation. The data processing adopts the method of coloring point clouds based on orthophotos. 2) Using the profile function, create a profile in the overlapping area of ​​each flight strip on the main interface and observe the profile view to see if there is obvious layering. If the point cloud is found to be layered, the point cloud calculation and flight strip adjustment operation need to be performed again. In addition, the quality inspection function in the intelligent laser module of the drone manager can be used to check the height difference of the point cloud between flight strips, which can also reflect whether the point cloud is layered. 3) Coordinate transformation: The required point cloud coordinate system is CGCS2000 National Geodetic Coordinate System, Gaussian 3-degree zone projection, central meridian 117°, and 1985 National Height Datum. Coordinate transformation is performed after point cloud adjustment and colorization, mainly involving parameter calculation and coordinate transformation. The residuals dN, dE, and dU of the seven parameters must be controlled within tolerance limits. Detailed steps are as follows: (1) The initial point cloud is in the WGS84 system. Now, to obtain the point cloud in the CGCS2000 national geodetic coordinate system, it is necessary to perform parameter transformation between the two coordinate systems. Using the three parameter-finding points measured in the field, and knowing the three pairs of coordinates of the three points in the WGS84 coordinate system and the CGCS2000 national geodetic coordinate system, parameter calculation is performed to obtain the seven parameters of the two coordinate systems. Check whether the residual information dN, dE and dU meet the requirements. (2) Configure the coordinate transformation parameters of the survey area, add the Bursa seven parameters obtained from the parameter calculation, complete the transformation parameter configuration, and transform the point cloud from the source coordinate system (WGS84 coordinate system) to the target coordinate system (CGCS2000 national geodetic coordinate system). 4) Standard point cloud data output: After the above preprocessing operations, UAV-borne LiDAR point cloud data can be generated for production. The resulting point cloud data is generally stored in "*.las" format. Ground point cloud data and image data can be used to generate 3D models and digital products such as DEM, DSM, and DOM as needed.

9. A method for precise monitoring of surface settlement during grouting in goaf areas according to claim 1, characterized in that, In step two, point cloud data processing follows the initial point cloud data processing. Further processing of the acquired point cloud data is required based on actual needs. The detailed steps are as follows: 1) Point cloud denoising Visualizing the standard format LAS point cloud of the mining area reveals a large number of noise points. Noise point clusters exist both above and below the surface point cloud. To ensure the accuracy of data products such as DEM generated from the point cloud data, it is necessary to remove the noise point clusters. The point cloud denoising function in the UAV Manager's Smart Point Cloud module is used to denoise the point cloud data of the entire project area. 2) Point cloud filtering and classification Point cloud filtering refers to the process of removing non-ground points (vegetation, buildings, bridges, power facilities, etc.) from point cloud data in order to extract digital elevation models. This process is known as LiDAR data filtering. Point cloud classification refers to the process of categorizing non-ground points for use in ground feature extraction and 3D reconstruction of buildings. 3) Digital product generation After using the drone management smart point cloud module to denoise and automatically classify standard point clouds, high-precision laser point cloud data suitable for digital product production is obtained. The data editing function in the smart point cloud module can generate a Digital Model (DSM) from the denoised point cloud, and further generate a Digital Image Model (DEM) from the automatically classified ground point cloud. For areas with poor automatic classification, manual classification can be performed through region smoothing, online classification, and offline classification. The resulting DSM and DEM digital products are shown. The DSM is generated from the denoised and classified point cloud. A DSM is generated from the obtained ground point cloud after point cloud denoising and automatic classification. 4) Data accuracy check, mainly including Point cloud data density and point cloud data plane and elevation accuracy checks are performed. The point density measurement function of the drone manager's intelligent laser module can measure the point density of the point cloud data. The side length of the measurement box can be set to measure the point density in different areas. Three types of areas were selected in the project area, and their point density was measured respectively. The three types of areas are dense building area, dense vegetation area, and surface area. The accuracy check function can check the elevation accuracy of the point cloud. By importing the check point file, it is compared with the point cloud data obtained by point cloud calculation to obtain the elevation difference between the point cloud data and the check point coordinates, thus realizing the accuracy check of the point cloud data.

10. A method for precise monitoring of surface settlement during grouting in goaf areas according to claim 1, characterized in that, Step five involves analyzing the coupling relationship between the residual deformation of the goaf and the compressive deformation of the building foundation. The detailed steps are as follows: 1) Calculation method for residual surface deformation under grouting conditions in goaf areas (1) Establish the theoretical model of probability integral method Grouting and filling of goaf areas can mitigate residual deformation of the surface. The surface movement is a relatively gentle residual settlement, which is basin-shaped with a larger center and smaller edges above the goaf. The morphology of the subsided basin still conforms to the normal curve distribution law. The mathematical model of probability integral method can be used to calculate the residual surface deformation. By establishing the theoretical model of probability integral method and through mathematical derivation, the expression of the subsidence value We(x,z) at point (x,z) caused by unit mining can be obtained. In the formula: W e (x,z)——Subsidence caused by unit mining, in mm; r z —Primarily affects the radius, in meters. For the surface, z equals the mining depth, and A is a constant representing the size of the medium unit, then r z Let r be a constant. z Let r be a constant, then the expression for a subsidence basin, a surface unit, is: (2) Calculation of residual surface deformation after grouting and filling of goaf area Applying the probability integral method formula, let the mining thickness at a point (s,t) within the goaf calculation area be Md(s,t). Then the subsidence at any point (x,y) on the surface is: In the formula: l A l B —The length and width of the goaf calculation area, in meters; M d — Coal seam mining thickness, in meters; q c —Residual surface subsidence coefficient of the goaf after grouting and filling; α—Coal seam dip angle; r — the main influence radius, in meters. Then, the tilt, curvature, horizontal displacement, and horizontal deformation at any point on the Earth's surface in the given direction φ are: (3) Analysis of residual surface settlement calculation parameters The residual surface deformation after grouting and filling of the goaf can be calculated by back-calculating the residual subsidence coefficient based on the goaf filling rate. (a) Coefficient of residual surface subsidence qc after grouting and filling of goaf area The residual surface subsidence coefficient is closely related to the subsidence coefficient during the duration of surface movement and the duration after mining ends. The calculation method for the residual surface movement deformation subsidence coefficient is as follows: q—Subsidence coefficient during the duration of surface movement; q1—Residual surface subsidence coefficient under normal mining conditions; k – Adjustment coefficient, typically ranging from 0.5 to 1.0; t — time remaining until the end of mining; The surface movement and deformation under grouting and filling conditions in goaf areas can be calculated by back-calculating the residual subsidence coefficient using the goaf filling rate. Therefore, the residual surface subsidence coefficient under grouting and filling conditions in goaf areas can be calculated using the following formula: In the formula: q c —Residual surface subsidence coefficient of the goaf after grouting and filling; ρ c —Ground filling rate; ρ y — Compaction rate of the filling material in the goaf. The compaction rate of the goaf filling body can be obtained through compression tests. The calculation method can be based on formula 11, and is generally taken as about 0.

1. In the formula: V – Volume of the filling material injected into the goaf, in meters (m³) 3 ; V' — Volume of the infill material after compression, in meters (m³) 3 (b) Inflection point offset coefficient S / H When the overlying rock types are hard, medium hard, and soft, the S / H values ​​when the surface is basically stable are 0.31–0.43, 0.08–0.30, and 0–0.07, respectively. The inflection point offset coefficient for calculating residual surface subsidence in goaf areas can be compared with the parameters when the surface is basically stable. (c) Main influencing angle tangent tanβ The tangent tanβ does not change significantly with the extension of time after mining ends. Its value is taken with reference to tanβ of conventional subsidence deformation. When the overlying rock type is hard, medium hard, and soft, the values ​​are 1.2~1.91, 1.92~2.40, and 2.41~3.54, respectively. (d) Horizontal movement coefficient b The horizontal movement coefficient b is mainly related to the properties of the overlying strata and has little to do with other factors. When calculating the horizontal movement coefficient b for residual surface settlement, the value of conventional settlement deformation can be referenced. The value of the horizontal movement coefficient ranges from 0.2 to 0.4, and is generally 0.

3. (e) Mining influence propagation angle θ The propagation angle θ of the mining impact depends on the coal seam dip angle and the mining impact propagation coefficient K. The mining impact propagation coefficient K is mainly related to the lithology of the overlying strata and the mining depth of the coal seam, while the impact of grouting and filling in the goaf is limited. To calculate the mining impact propagation angle θ during residual surface settlement, the value of θ for conventional subsidence deformation can be referenced: θ = 90° - kα, where α is the coal seam dip angle, and k is taken as 0.7–0.8, 0.6–0.7, and 0.5–0.6 for hard, medium-hard, and weak overlying strata, respectively. 2) Calculation model for soil compression settlement under building loads When constructing engineering projects above a goaf, the newly built buildings and structures are affected not only by the residual deformation of the goaf but also by the compression and settlement deformation of the foundation soil under the load of the building itself. Therefore, the residual surface deformation above the goaf should include the superposition of the residual movement deformation of the goaf and the compression and settlement of the soil under the load. (1) Spatial stress analysis of foundation soil under load The influence range of a rectangular load (length and width a×b) on the underlying foundation soil is approximately an ellipsoid. For ease of calculation, the maximum influence range of the soil can be considered as a cuboid with length, width, and height la, lb, and Hz, where la and lb are the length and width of the soil affected by the load, and Hz is the depth of the load influence. Therefore, the line connecting the critical points where the load influence depth is 0 along the center lines of the long and short sides of the foundation is taken as the lengths of la and lb, and Hz is taken as the maximum value of the load influence depth. A point M with coordinates (x, y, z) is taken within the affected space of the soil. The rectangular foundation is subjected to a uniformly distributed load P. A unit cell is taken at point M. The additional vertical and horizontal stresses at the unit cell can be calculated according to Equation 12. In the formula: σ fz(x,y,z) σ fx(x,y,z) σ fy(x,y,z) —Additional stress in the x, y, and z directions at point (x, y, z); σ z(x,y,z) σ x(x,y,z) σ y(x,y,z) —Additional stress coefficients in the x, y, and z directions at point (x, y, z) (2) Layered settlement model of foundation soil The foundation soil exhibits the characteristics of natural soil layers, and considering the compressibility of the soil, the layered summation method from soil mechanics is applied to calculate the foundation settlement. This model makes the following assumptions: (a) Each layer in the foundation undergoes vertical compressive deformation without lateral expansion. (b) The layer thickness should not be too large, generally 1 to 2 m. The stress of each layer from top to bottom varies. When calculating, the average value of the stress at the top and bottom of the layer is taken approximately. (3) Calculate the soil deformation within the load influence range. The settlement at a point (x, y, 0) on the surface of the foundation soil is the cumulative compressive settlement of all soil elements along the straight line from that point to the depth z = Hz. The settlement value of the element at point M(x, y, z) can be calculated using Equation 13. (X,Y,Z) Integrating over z = [0, Hz], the settlement at point (x, y, 0) can be calculated using Equation 14: In the formula: dS (X,Y,Z) —Vertical settlement of the soil element at point (x,y,z); E S(X,Y,Z) —The compression modulus of the soil element at point (x,y,z); S (x,y,0) — Settlement at point (x,y) on the surface of the foundation soil (z=0). Because the foundation soil has the characteristic of natural soil layering, it can be divided into m1 soil layers from the z=0 plane to the Z=Hz plane, then S (x,y,0) Formula 14 can also be used for calculation: In the formula: —The average additional stress in the i-th soil layer below point (x,y,0); —The average additional stress coefficient of the i-th soil layer below point (x,y,0); Δz — the thickness of the i-th soil layer below point (x,y,0); E s —The compressibility modulus of the i-th soil layer below point (x,y,0); m1—Number of soil layers within the depth range affected by the load; p—mean pressure at the base, in kg / m². The settlement of building foundations under load can be considered as the superposition of the compression settlement of the foundation soil and the residual surface settlement, i.e., S 总 =S d +S w When calculating the settlement of the two, the origins of the selected coordinate systems are different. By transforming the coordinate system for calculating the compression of the foundation soil to the coordinate system for calculating the residual deformation of the goaf, the synergistic relationship between the residual settlement of the goaf and the compression deformation of the building foundation can be calculated according to Formula 15. The settlement and horizontal movement of the foundation soil are mainly due to the coordinated settlement of residual deformation in the goaf and compressive deformation of the foundation soil. Therefore, the tilt, horizontal movement, curvature, and horizontal deformation in the given direction p at any point on the surface of the foundation soil can be expressed by the following formula. (4) Calculation of settlement and deformation of building foundation under grouting and filling conditions in goaf area The settlement of building foundations under grouting and filling conditions in goaf areas mainly includes two parts: residual surface settlement in the goaf area and compressive settlement of the foundation soil. Specific calculations are given in formulas 17-20. In the calculation, the residual surface settlement coefficient can be calculated using formula 9. The compressive settlement of the foundation soil is calculated using the layered summation method, thus allowing the calculation of the settlement of 1 / 4 of the rectangular space of the foundation soil. The settlement at a point on the surface of the foundation soil is the cumulative compressive settlement of all soil units along the straight line from that point to the depth of load influence. Layering is required according to the soil type and layer thickness, and the layer thickness should not be too large (not exceeding 3m). Formula 14 is used to calculate the settlement at any point on the soil surface.

Images