Underground mining settlement monitoring method and system based on InSAR technology

Abstract

The invention relates to the technical field of surveying and mapping and remote sensing, and discloses an underground mining settlement monitoring method and system based on an InSAR technology. The method comprises the following steps: acquiring a multi-temporal SAR image, GNSS observation data, a DEM and a mining log; generating a differential interferogram sequence; constructing an atmospheric phase screen model by using a GNSS for correction; establishing a theoretical settlement prior model based on the mining parameters; establishing a theoretical settlement priori model by using the mining parameters as a constraint to embed an improved minimum cost flow phase unwrapping algorithm; and performing space-time filtering and time sequence stacking on the unwrapping phase, and generating a ground surface three-dimensional deformation field time sequence in combination with singular value decomposition and adaptive deformation model fitting. The system comprises a multi-source data acquisition module, a differential interferogram generation module, an atmospheric correction module, a prior modeling module, a constraint unwrapping module, a deformation reconstruction module and the like. Millimeter-level precision, high-frequency and large-range mining area settlement dynamic monitoring is realized through multi-source fusion and physical mechanism constraint.

Description

Technical Field

[0001] This invention belongs to the field of surveying and remote sensing technology, specifically relating to a method and system for monitoring underground mining subsidence based on InSAR technology. Background Technology

[0002] With the rapid advancement of urbanization and resource development, surface subsidence caused by underground mining activities has become a significant geological hazard threatening the ecological security and infrastructure stability of mining areas. Synthetic Aperture Radar Interferometry (InSAR) technology, with its advantages of wide coverage, high precision, and non-contact operation, has been widely applied in the field of surface deformation monitoring, capable of acquiring millimeter-level temporal information on surface subsidence. However, InSAR can only reflect the macroscopic results of surface deformation and cannot reveal the underlying hydrogeological driving mechanisms of subsidence. It struggles to distinguish subsidence responses caused by different factors such as goaf collapse, aquifer depletion, or tectonic activity, thus limiting its in-depth application in risk warning and causal tracing.

[0003] Dynamic changes in groundwater levels are a significant factor inducing or exacerbating land subsidence. Drainage of aquifers increases effective stress, leading to soil compression and surface subsidence. Traditional monitoring methods rely heavily on sparsely deployed groundwater level observation wells, which, while providing localized water level data, have limited spatial coverage, low update frequency, and mismatched spatiotemporal resolution with InSAR deformation fields, making it difficult to establish a dynamic coupling relationship between the two. Existing technologies typically involve simple overlaying or post-hoc correlation analysis of InSAR deformation data and hydrological data, lacking deep integration in terms of temporal synchronization, spatial correspondence, and consistency of physical mechanisms, thus failing to construct predictive subsidence-water level response models.

[0004] Currently, no existing technology can achieve accurate spatiotemporal registration of InSAR surface deformation sequences and automatic high-frequency groundwater level monitoring data, nor can it quantify the contribution of aquifer dewatering to land subsidence, or set dynamic thresholds based on real-time water level changes to trigger subsidence risk warnings. This technological gap makes it difficult for mine managers to identify high-risk areas before subsidence occurs, and also hinders the implementation of differentiated prevention and control strategies for different hydrogeological conditions, thus restricting the initiative and timeliness of ecological restoration and infrastructure protection. There is a need for an underground mining subsidence monitoring method and system that integrates multi-source heterogeneous monitoring data and possesses mechanistic interpretation capabilities and early warning functions. Summary of the Invention

[0005] This invention provides a method and system for monitoring subsidence in underground mining based on synthetic aperture radar interferometry, aiming to solve the problems of insufficient subsidence inversion accuracy, long monitoring cycle, low spatial resolution, and weak dynamic response capability in existing technologies due to atmospheric delay errors, spatiotemporal decoherence effects, and single data sources. This invention achieves high-precision, high-frequency, and large-scale continuous monitoring of millimeter-level subsidence fields in mining areas by constructing a multi-source heterogeneous remote sensing data fusion framework, introducing an external physical constraint model, and designing adaptive phase unwrapping and deformation temporal reconstruction algorithms.

[0006] This invention provides a method for monitoring subsidence in underground mining based on InSAR technology, comprising: Acquire a dataset of multi-temporal synthetic aperture radar images covering the target mining area; Simultaneously acquire global navigation satellite system observation data, digital elevation model data, and mining progress logs that are time-matched with the multi-temporal synthetic aperture radar image dataset; Accurate registration and baseline estimation are performed on the multi-temporal synthetic aperture radar image dataset to generate a differential interferogram sequence; An atmospheric phase screen model is constructed using observation data from the Global Navigation Satellite System, and atmospheric delay phase correction is performed on the differential interferogram sequence. Based on the digital elevation model data and the mining progress log, a theoretical subsidence prior model caused by mining is established. The theoretical settlement prior model is used as an external constraint and embedded into the improved minimum cost flow phase unwrapping algorithm to generate the unwrapped interference phase sequence. The unwrapped interference phase sequence is subjected to spatiotemporal filtering and time-series stacking processing. The main deformation mode is extracted by singular value decomposition and parameter fitting is performed by combining linear and nonlinear deformation rate models to finally generate a three-dimensional deformation field time series product of the Earth's surface.

[0007] Preferably, synchronously acquiring Global Navigation Satellite System (GNSS) observation data that is time-matched to the multi-temporal synthetic aperture radar (SAR) image dataset includes: No fewer than three permanent global navigation satellite system monitoring stations will be set up in the mining area to record the carrier phase and pseudorange observations within two hours before and after the imaging time of each synthetic aperture radar image. The three-dimensional coordinates and covariance matrix of each monitoring station at the imaging time will be calculated by a precise single-point positioning algorithm.

[0008] Preferably, constructing an atmospheric phase screen model using observation data from the Global Navigation Satellite System includes: The total zenith delay calculated by the global navigation satellite system is converted into atmospheric path delay phase along the radar line of sight. Using the Kriging interpolation method, with the locations of the Global Navigation Satellite System monitoring stations as control points, a two-dimensional spatial distribution field of atmospheric phase covering the entire differential interferogram region is generated. The atmospheric phase two-dimensional spatial distribution field is subtracted pixel by pixel from the original differential interferogram to complete the atmospheric phase correction.

[0009] Preferably, a priori model for theoretical settlement caused by mining is established, including: Based on the working face advancement position, mining thickness, mining depth and rock stratum movement angle parameters in the mining progress log, the expected subsidence value, tilt value and curvature value of any point on the surface are calculated using the probability integral method. The calculation results are rasterized into a priori settlement rate map with the same resolution as the synthetic aperture radar image, and a spatial weighting coefficient is assigned to it, which decreases exponentially with the distance from the center of the working face.

[0010] Preferably, the theoretical settlement prior model is used as an external constraint and embedded into the improved minimum cost flow phase unwrapping algorithm, including: When constructing the phase gradient network graph, the phase gradient corresponding to the theoretical settlement prior graph is used as the initial flow reference value of the edge; A prior constraint penalty term is added to the minimum cost flow optimization objective function, and its weight is dynamically adjusted by the local coherence coefficient.

[0011] Preferably, spatiotemporal filtering and time-series stacking processing are performed on the unwrapped interference phase sequence, including: A wavelet thresholding denoising method is used to perform spatial domain filtering on each unwrapped phase image; In the time dimension, a sliding window mid-range filter is applied to the phase sequence of the same geographic coordinates, with a window length of 5 images; Then, all the filtered phase maps are stacked in chronological order to form a three-dimensional phase cube.

[0012] Preferably, the singular value decomposition method is used to extract the main deformation mode, and the parameters are fitted by combining linear and nonlinear deformation rate models, including: Singular value decomposition is performed on the three-dimensional phase cube, and the top principal components with a cumulative contribution rate of more than 95% are retained. The principal component time series were fitted with the linear function, the quadratic polynomial function and the exponential decay function respectively using least squares fitting. The optimal deformation model type is automatically selected based on the principle of minimizing the root mean square of the fitting residuals. Finally, the parameters of the selected model are inverted into the deformation rate fields of the Earth's surface in the east, north, and vertical directions.

[0013] Preferably, during the generation of the differential interferogram sequence, after performing conjugate multiplication on each pair of interferometric combinations, the terrain phase component is simulated and removed using the digital elevation model data to obtain a differential interferogram containing only deformation, atmospheric, and noise information.

[0014] Preferably, the spatial weighting coefficient of the settlement rate prior map is defined as an exponential decay function with the current working face center as the origin, used to characterize the spatial distribution characteristics of the reliability of the theoretical model prediction.

[0015] This invention also provides an underground mining settlement monitoring system based on InSAR technology. The system uses the aforementioned InSAR-based underground mining settlement monitoring method to monitor underground mining settlement. The system includes: The multi-source remote sensing data acquisition unit is used to acquire multi-temporal synthetic aperture radar image datasets covering the target mining area, global navigation satellite system observation data, digital elevation model data, and mining progress logs of the mining area. The differential interferogram generation unit is used to perform accurate registration and baseline estimation on the multi-temporal synthetic aperture radar image dataset to generate a differential interferogram sequence. The atmospheric phase correction unit is used to construct an atmospheric phase screen model using the observation data of the Global Navigation Satellite System and to perform atmospheric delay phase correction on the differential interferogram sequence. The theoretical settlement prior modeling unit is used to establish a theoretical settlement prior model caused by mining based on the digital elevation model data and the mining progress log of the mining area; The constrained phase unwrapping unit is used to embed the theoretical settlement prior model as an external constraint into the improved minimum cost flow phase unwrapping algorithm to generate the unwrapped interference phase sequence. The deformation time series reconstruction unit is used to perform spatiotemporal filtering and time series stacking processing on the unwrapped interference phase sequence, extract the main deformation mode using the singular value decomposition method, and perform parameter fitting by combining linear and nonlinear deformation rate models, and finally generate a three-dimensional deformation field time series product of the earth's surface.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This invention constructs a new paradigm for settlement monitoring that combines physical mechanism-driven and data-driven approaches by integrating multi-orbit synthetic aperture radar imagery, atmospheric observation data from global navigation satellite systems, and mining engineering parameters.

[0017] 2. The atmospheric delay phase is accurately corrected using data from the Global Navigation Satellite System, eliminating spurious deformation signals caused by atmospheric disturbances in traditional methods; a theoretical settlement model based on the probability integral method is introduced as a priori constraint into the phase unwrapping process, improving the reliability and accuracy of unwrapping in low-coherence regions; and a complete characterization of the entire mining settlement process is achieved through singular value decomposition and adaptive selection mechanism of time series models.

[0018] 3. The three-dimensional deformation field time series product generated by this invention has optimized vertical monitoring accuracy, high horizontal resolution, and shortened time update cycle. It can provide timely and accurate decision support data for mine safety production, geological disaster early warning and ecological restoration. It overcomes the problems of existing technologies that rely on a single data source, ignore mining physical mechanisms, accumulate a lot of untangling errors and cannot distinguish between real deformation and noise interference. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the overall technical solution architecture of the present invention; Figure 2 This is a schematic diagram of the core principle framework of constrained phase untangling in this invention, which integrates the theoretical settlement prior model and the improved minimum cost flow algorithm. Figure 3 This is a logical flowchart of the multi-source heterogeneous remote sensing data acquisition and differential interferogram generation stage in this invention. Figure 4 This is a logical flowchart of the atmospheric phase correction and spatiotemporal filtering processing stage in this invention; Figure 5 This is a schematic diagram of the multi-level interaction relationship and data flow between the terminal system module and the external data source in this invention. Detailed Implementation

[0020] refer to Figures 1 to 5 This invention provides a method and system for monitoring subsidence in underground mining based on synthetic aperture radar interferometry. This method integrates multi-source heterogeneous remote sensing data with mining area engineering parameters to construct a high-precision deformation inversion framework constrained by physical mechanisms, enabling continuous, dynamic, and three-dimensional monitoring of the millimeter-level subsidence field on the mining area surface. The specific implementation of this invention will be described in detail below according to the steps explicitly defined in the invention description.

[0021] Step S1: Acquire a multi-temporal synthetic aperture radar (SAR) image dataset covering the target mining area. This step requires acquiring at least 15 single-look complex images from at least two SAR satellite platforms with different orbital orientations. The time interval between the selected images must be strictly controlled within 12 days to ensure effective capture of rapid deformation processes; the vertical baseline length must be maintained within 200 meters to suppress phase noise caused by terrain; the image coverage area must extend at least 5 kilometers beyond the actual boundary of the mining area to provide sufficient buffer space for subsequent atmospheric correction and edge effect processing. All acquired single-look complex images contain complete amplitude and phase information and record their precise imaging time, satellite position, orbital parameters, wavelength, and other metadata as the basic input for subsequent processing.

[0022] Step S2: Synchronously acquire GNSS observation data, digital elevation model (DEM) data, and mining progress logs that are time-matched to the multi-temporal synthetic aperture radar (SAR) image dataset. GNSS observation data acquisition relies on at least three permanent GNSS monitoring stations deployed within and around the target mining area. These stations continuously record the carrier phase and pseudorange raw observations for two hours before and after each SAR image imaging moment at a sampling frequency of at least once per second. Subsequently, these raw observations are processed using a precise single-point positioning algorithm to calculate the precise three-dimensional coordinates and covariance matrix of each monitoring station corresponding to each SAR image imaging moment. The DEM data uses a resolution of 10 meters or higher to accurately describe the surface topography. The mining progress log is provided by the mine production management system and records detailed engineering parameters such as the daily working face advancement position, mining thickness, mining depth, and rock strata movement angle.

[0023] Step S3: Accurate registration and baseline estimation are performed on the multi-temporal synthetic aperture radar (SAR) image dataset to generate a differential interferogram (DII) sequence. This process first resamples all single-view complex images to the same georeferenced grid, achieving sub-pixel accuracy. Then, based on precise satellite orbit parameters, the spatiotemporal baseline of each interferometric pair is calculated, and qualified interferometric pairs that meet the temporal and spatial decoherence thresholds are selected. For each grid-connected interferometric pair, the primary and secondary images are conjugate-multiplied to generate an original interferogram containing topographic phase, deformation phase, atmospheric phase, and noise phase. Next, high-precision digital elevation model (DEM) data is used to simulate and remove the topographic phase component, ultimately obtaining a DPI sequence containing only deformation, atmospheric, and noise information.

[0024] Step S4 involves constructing an atmospheric phase screen model using the observation data from the Global Navigation Satellite System (GNSS) and performing atmospheric delay phase correction on the differential interferogram (DIC) sequence. This step converts the total zenith delay calculated by the GNSS into atmospheric path delay phase along the line-of-sight direction of the synthetic aperture radar (SAR) using a mapping function. Due to the limited number of GNSS monitoring stations, Kriging interpolation is employed, using their three-dimensional coordinates and corresponding atmospheric phase values ​​as control points to generate a continuous and smooth two-dimensional spatial distribution field of atmospheric phase within the entire DIC coverage area. This distribution field constitutes the atmospheric phase screen model. The correction process is achieved by subtracting this atmospheric phase screen model pixel-by-pixel from the original DIC, thereby eliminating atmospheric disturbance phase caused by spatiotemporal variations in water vapor content and improving the signal-to-noise ratio of the true deformation signal in the DIC.

[0025] Step S5: Based on the digital elevation model data and the mining progress log, a theoretical subsidence prior model caused by mining is established. This model is constructed using the classic probabilistic integral method. First, based on the working face geometric parameters (including advance position, mining thickness, and mining depth) and rock mechanics parameters (mainly the rock strata movement angle) in the mining progress log, the expected subsidence, tilt, and curvature values ​​of any point on the surface under the current mining conditions are calculated. This calculation process fully considers the spatial distribution characteristics of mining activities and the physical laws of rock strata movement. Subsequently, the calculated continuous subsidence field is rasterized to make its spatial resolution completely consistent with the synthetic aperture radar image, forming a subsidence rate prior map. To reflect the reliability of prior information changing with spatial location, a spatial weight coefficient is assigned to this prior map. This spatial weight coefficient is defined as an exponential decay function with the current working face center as the origin. The closer to the working face, the higher the weight, indicating that the theoretical model's prediction in this area is more reliable; conversely, in stable areas far from the working face, the weight approaches 0 to avoid introducing unnecessary bias.

[0026] Step S6: The theoretical settlement prior model is used as an external constraint and embedded into the improved minimum cost flow phase unwrapping algorithm to generate the unwrapped interference phase sequence. Traditional phase unwrapping algorithms are prone to errors in low coherence regions, leading to distortion of deformation information. This invention guides the unwrapping process by introducing physical prior knowledge. Specifically, during the construction of the phase gradient network graph, the phase gradient value corresponding to the theoretical settlement prior graph is used as the initial flow reference value for each edge in the network. More importantly, a prior constraint penalty term is added to the optimization objective function of the minimum cost flow algorithm. The mathematical expression of this prior constraint penalty term is:

[0027] The untangling phase to be solved. The phase is derived from the theoretical settlement prior model. Represents the gradient operator, These are dynamic weighting coefficients, whose values ​​are determined by the local coherence coefficient. The decision is made when the local coherence is higher than 0.7, indicating that the interference signal quality in the low-coherence region is high and the prior information is highly trusted. The value is 0.8; when the local coherence is below 0.4, it indicates that the signal in that low-coherence region is severely contaminated by noise and mainly depends on the data itself. The value is 0.2; it is between 0.4 and 0.7. The method is determined by linear interpolation. Through this adaptive mechanism, the algorithm can fully utilize the accuracy of the prior model in the high coherence region and avoid being misled by unreliable priors in the low coherence region, thereby globally improving the robustness and accuracy of the untangling results.

[0028] Step S7: The unwrapped interferometric phase sequence is subjected to spatiotemporal filtering and temporal stacking. Singular value decomposition (SVD) is used to extract the principal deformation modes, and linear and nonlinear deformation rate models are combined for parameter fitting to ultimately generate a three-dimensional deformation field time series product. This step first performs spatial domain filtering on each unwrapped phase image, using wavelet thresholding to suppress residual high-frequency noise. Then, in the temporal dimension, for each fixed geographic coordinate point, a sliding window midpoint filter is applied to its phase time series, with a window length set to 5 images to further smooth the time series and remove isolated outliers. After the above dual filtering, all phase images are stacked in chronological order to form a three-dimensional phase cube, with dimensions of longitude, latitude, and time. Singular value decomposition is performed on this three-dimensional phase cube to obtain a set of mutually orthogonal spatiotemporal principal components. The top few principal components with a cumulative contribution rate exceeding 95% are retained; these principal components represent the main spatiotemporal modes of regional deformation. Next, the time series of each principal component is fitted with least-squares models of three candidate deformation models: linear functions, quadratic polynomial functions, and exponentially decaying functions. By comparing the root mean square of the fitting residuals of each model, the optimal deformation model type is automatically selected. Finally, the parameters of the selected optimal model are inverted and decomposed into three geographical directions: east, north, and vertical, generating a complete and high-precision three-dimensional deformation field time series product.

[0029] Corresponding to the above method, the present invention also provides an underground mining settlement monitoring system based on synthetic aperture radar interferometry. This underground mining settlement monitoring system, as an integrated data processing and analysis platform, has an internal structure that strictly corresponds to the steps of the aforementioned method.

[0030] The system includes a multi-source remote sensing data acquisition unit. This unit is responsible for coordinating and scheduling external data sources. Specifically, it is configured to acquire no fewer than 15 single-view complex images from at least two synthetic aperture radar (SAR) satellite platforms with different orbital orientations, ensuring that the time intervals, vertical baselines, and coverage meet the aforementioned requirements. Simultaneously, this multi-source remote sensing data acquisition unit manages no fewer than three permanent GNSS monitoring stations deployed within the mining area, controlling them to acquire raw observation data at a frequency of no less than once per second, and receiving digital elevation model data and mining progress logs from the mine production management system.

[0031] The system includes a differential interferogram generation unit. This unit receives a multi-temporal synthetic aperture radar image dataset from the multi-source remote sensing data acquisition unit, performs precise image registration, interferometric pair selection, conjugate multiplication, and terrain phase removal operations, and finally outputs a sequence of differential interferograms.

[0032] The atmospheric phase correction unit is connected to the differential interferogram generation unit. This atmospheric phase correction unit and differential interferogram generation unit receive observation data from the Global Navigation Satellite System and differential interferogram sequences, perform the conversion from zenith total delay to line-of-sight phase, generate an atmospheric phase screen model through Kriging interpolation, and complete the atmospheric phase correction of the differential interferogram, outputting the corrected differential interferogram.

[0033] The theoretical settlement prior modeling unit receives digital elevation model data and mining progress logs. This unit has a built-in probability integral method calculation engine that can calculate the theoretical surface settlement field in real time based on the input engineering parameters, perform rasterization and spatial weight assignment, and output a prior map of settlement rates.

[0034] The constrained phase unwrapping unit receives the differential interferogram sequence output by the atmospheric phase correction unit and the prior map output by the theoretical settlement prior modeling unit. This constrained phase unwrapping unit implements the aforementioned improved minimum cost flow phase unwrapping algorithm, and can dynamically adjust the prior constraint weights to generate a high-precision unwrapped interferometric phase sequence.

[0035] The deformation time-series reconstruction unit, as the final output module of the system, receives the unwrapped interference phase sequence. This unit integrates a series of functions, including wavelet denoising, sliding window midpoint filtering, three-dimensional phase cube construction, singular value decomposition, multi-model fitting, and optimal model selection. Its final output is a time-series product of the deformation rate field of the Earth's surface in the east, north, and vertical directions, which can be directly used for downstream applications such as mine safety assessment, geological disaster early warning, and ecological restoration planning.

[0036] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0037] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for monitoring underground mining subsidence based on InSAR technology, characterized in that, include: Acquire a dataset of multi-temporal synthetic aperture radar images covering the target mining area; Simultaneously acquire global navigation satellite system observation data, digital elevation model data, and mining progress logs that are time-matched with the multi-temporal synthetic aperture radar image dataset; Accurate registration and baseline estimation are performed on the multi-temporal synthetic aperture radar image dataset to generate a differential interferogram sequence; An atmospheric phase screen model is constructed using observation data from the Global Navigation Satellite System, and atmospheric delay phase correction is performed on the differential interferogram sequence. Based on the digital elevation model data and the mining progress log, a theoretical subsidence prior model caused by mining is established. The theoretical settlement prior model is used as an external constraint and embedded into the improved minimum cost flow phase unwrapping algorithm to generate the unwrapped interference phase sequence. The unwrapped interference phase sequence is subjected to spatiotemporal filtering and time-series stacking processing. The main deformation mode is extracted by singular value decomposition and parameter fitting is performed by combining linear and nonlinear deformation rate models to finally generate a three-dimensional deformation field time series product of the Earth's surface.

2. The method for monitoring underground mining subsidence based on InSAR technology according to claim 1, characterized in that, Synchronously acquire Global Navigation Satellite System (GNSS) observation data that is time-matched to the multi-temporal synthetic aperture radar (SAR) image dataset, including: No fewer than three permanent global navigation satellite system monitoring stations will be set up in the mining area to record the carrier phase and pseudorange observations within two hours before and after the imaging time of each synthetic aperture radar image. The three-dimensional coordinates and covariance matrix of each monitoring station at the imaging time will be calculated by a precise single-point positioning algorithm.

3. The method for monitoring underground mining subsidence based on InSAR technology according to claim 2, characterized in that, Constructing an atmospheric phase screen model using observation data from the aforementioned Global Navigation Satellite System includes: The total zenith delay calculated by the global navigation satellite system is converted into atmospheric path delay phase along the radar line of sight. Using the Kriging interpolation method, with the locations of the Global Navigation Satellite System monitoring stations as control points, a two-dimensional spatial distribution field of atmospheric phase covering the entire differential interferogram region is generated. The atmospheric phase two-dimensional spatial distribution field is subtracted pixel by pixel from the original differential interferogram to complete the atmospheric phase correction.

4. The method for monitoring underground mining subsidence based on InSAR technology according to claim 3, characterized in that, Establish a priori model for theoretical settlement caused by mining, including: Based on the working face advancement position, mining thickness, mining depth and rock stratum movement angle parameters in the mining progress log, the expected subsidence value, tilt value and curvature value of any point on the surface are calculated using the probability integral method. The calculation results are rasterized into a priori settlement rate map with the same resolution as the synthetic aperture radar image, and a spatial weighting coefficient is assigned to it, which decreases exponentially with the distance from the center of the working face.

5. The method for monitoring underground mining subsidence based on InSAR technology according to claim 4, characterized in that, The theoretical settlement prior model is used as an external constraint and embedded into the improved minimum cost flow phase unwrapping algorithm, including: When constructing the phase gradient network graph, the phase gradient corresponding to the theoretical settlement prior graph is used as the initial flow reference value of the edge; A prior constraint penalty term is added to the minimum cost flow optimization objective function, and its weight is dynamically adjusted by the local coherence coefficient.

6. The method for monitoring underground mining subsidence based on InSAR technology according to claim 5, characterized in that, The unwrapped interference phase sequence is subjected to spatiotemporal filtering and time-series stacking processing, including: A wavelet thresholding denoising method is used to perform spatial domain filtering on each unwrapped phase image; In the time dimension, a sliding window mid-range filter is applied to the phase sequence of the same geographic coordinates, with a window length of 5 images; Then, all the filtered phase maps are stacked in chronological order to form a three-dimensional phase cube.

7. The method for monitoring underground mining subsidence based on InSAR technology according to claim 6, characterized in that, The singular value decomposition method is used to extract the principal deformation modes, and parameter fitting is performed by combining linear and nonlinear deformation rate models, including: Singular value decomposition is performed on the three-dimensional phase cube, and the top principal components with a cumulative contribution rate of more than 95% are retained. The principal component time series were fitted with the linear function, the quadratic polynomial function and the exponential decay function respectively using least squares fitting. The optimal deformation model type is automatically selected based on the principle of minimizing the root mean square of the fitting residuals. Finally, the parameters of the selected model are inverted into the deformation rate fields of the Earth's surface in the east, north, and vertical directions.

8. The method for monitoring underground mining subsidence based on InSAR technology according to claim 7, characterized in that, In the process of generating the differential interferogram sequence, after performing conjugate multiplication on each pair of interferometric combinations, the terrain phase component is simulated and removed using the digital elevation model data to obtain a differential interferogram that contains only deformation, atmospheric and noise information.

9. The method for monitoring underground mining subsidence based on InSAR technology according to claim 8, characterized in that, The spatial weighting coefficient of the settlement rate prior map is defined as an exponential decay function with the current working face center as the origin, used to characterize the spatial distribution characteristics of the reliability of the theoretical model prediction.

10. A subsidence monitoring system for underground mining based on InSAR technology, characterized in that, The method for monitoring underground mining settlement based on InSAR technology, as described in any one of claims 1 to 9, is used to implement underground mining settlement monitoring. The system includes: The multi-source remote sensing data acquisition unit is used to acquire multi-temporal synthetic aperture radar image datasets covering the target mining area, global navigation satellite system observation data, digital elevation model data, and mining progress logs of the mining area. The differential interferogram generation unit is used to perform accurate registration and baseline estimation on the multi-temporal synthetic aperture radar image dataset to generate a differential interferogram sequence. The atmospheric phase correction unit is used to construct an atmospheric phase screen model using the observation data of the Global Navigation Satellite System and to perform atmospheric delay phase correction on the differential interferogram sequence. The theoretical settlement prior modeling unit is used to establish a theoretical settlement prior model caused by mining based on the digital elevation model data and the mining progress log of the mining area; The constrained phase unwrapping unit is used to embed the theoretical settlement prior model as an external constraint into the improved minimum cost flow phase unwrapping algorithm to generate the unwrapped interference phase sequence. The deformation time series reconstruction unit is used to perform spatiotemporal filtering and time series stacking processing on the unwrapped interference phase sequence, extract the main deformation mode using the singular value decomposition method, and perform parameter fitting by combining linear and nonlinear deformation rate models, and finally generate a three-dimensional deformation field time series product of the earth's surface.

Images