Space-air-ground integrated mine area geological deformation millimeter level monitoring method and system

Abstract

This invention discloses a millimeter-level monitoring method and system for geological deformation in mining areas based on an integrated space-air-ground system, relating to the field of geological monitoring technology. The method includes: acquiring time-series synthetic aperture radar (SAR) images, GNSS observation data, and subsurface microseismic data; inverting the physical property parameters of the surface cover medium based on GNSS observation data, calculating the medium phase delay, and compensating for the original interferometric encapsulation phase to obtain the true geometric encapsulation phase; performing spatial clustering based on subsurface microseismic data to extract fault rupture surfaces and generate discontinuous masks, which are then superimposed onto the true geometric encapsulation phase to block cross-fault phase unwrapping paths, dividing the mining area into multiple topologically isolated blocks and performing independent phase unwrapping on each; acquiring the absolute deformation projection values ​​of GNSS stations within each block, combining the unwrapped phase to calculate the reference intercept, and reconstructing the continuous absolute deformation field of the mining area; this invention achieves high-precision monitoring of mining area deformation through medium phase compensation and fault topological constraints.

Description

Technical Field

[0001] This invention relates to the field of geological monitoring technology, and in particular to a method and system for millimeter-level monitoring of geological deformation in mining areas based on an integrated air-space-ground system. Background Technology

[0002] Monitoring geological deformation in mining areas is a key technology for ensuring safe production and environmental management in mines. It primarily involves using integrated space-air-ground methods, such as Synthetic Aperture Radar Interferometry (InSAR), Global Navigation Satellite System (GNSS), and underground microseismic monitoring, to acquire displacement information at the surface and deep layers of the mining area. This technology is widely applied in scenarios such as coal mine subsidence area management, early warning of hidden fault activity, and slope stability evaluation, aiming to achieve high-precision perception of the complex geological evolution processes in mining areas.

[0003] Existing monitoring approaches are typically based on InSAR phase unwrapping algorithms, assuming that surface deformation satisfies phase continuity conditions in space and ignoring the influence of fine-grained media environments beyond the atmosphere. However, in specific mining areas with wind-blown sand cover and activated hidden faults, the dielectric properties of the sand layer vary with thickness and dielectric constant, resulting in additional propagation delays. Furthermore, severe fault displacement can cause discontinuous steps in surface deformation. The unwrapping path planning characteristics based on the continuous phase assumption in existing technologies fail to effectively remove spurious deformation phases caused by wind-blown sand, and are prone to accumulating unwrapping errors in cross-fault regions, severely impacting the absolute accuracy of deformation monitoring.

[0004] Therefore, a method and system for monitoring geological deformation in mining areas at the millimeter level based on an integrated air-space-ground system is proposed. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method and system for millimeter-level monitoring of geological deformation in mining areas based on an integrated air-space-ground system.

[0006] To achieve the above objectives, the technical solution of the present invention is as follows:

[0007] According to one aspect of this application, a millimeter-level monitoring method for geological deformation in mining areas based on an integrated air-space-ground system is provided, comprising: Acquire multi-source observation data of the target mining area, including time-series synthetic aperture radar imagery, global navigation satellite system observation data, and underground microseismic data; The physical property parameters of the surface cover medium of the target mining area are inverted based on observation data from the Global Navigation Satellite System. The phase delay of the medium is calculated based on the physical property parameters and the radar wave propagation mechanism. The original interferometric wrapping phase of the time-series synthetic aperture radar image is compensated to obtain the true geometric wrapping phase. Spatial clustering based on underground microseismic data is used to extract fault rupture surfaces and generate discontinuous masks. The discontinuous masks are superimposed on the real geometrically wrapped phase to block the phase unwrapping integration path across faults. The target mining area is divided into multiple topologically isolated blocks, and each block is independently unwrapped. Obtain the absolute deformation projection values ​​of the global navigation satellite system stations within each block, calculate the reference intercept of the corresponding block by combining the unwrapped phase of each block, and reconstruct the continuous absolute deformation field of the target mining area based on the reference intercept.

[0008] According to another aspect of this application, a millimeter-level monitoring system for geological deformation in mining areas based on an integrated air-space-ground system is provided, comprising: The data acquisition module is used to acquire multi-source observation data of the target mining area, including time-series synthetic aperture radar imagery, global navigation satellite system observation data, and underground microseismic data. The phase compensation module is used to invert the physical property parameters of the surface covering medium of the target mining area based on observation data from the Global Navigation Satellite System, and calculate the phase delay of the medium based on the physical property parameters and the radar wave propagation mechanism. It compensates the original interferometric wrapping phase of the time-series synthetic aperture radar image to obtain the true geometric wrapping phase. The phase unwrapping module is used to perform spatial clustering based on underground microseismic data to extract fault rupture surfaces and generate discontinuous masks. The discontinuous masks are superimposed on the real geometric wrapping phase to block the phase unwrapping integration path across faults, dividing the target mining area into multiple topologically isolated blocks, and performing independent phase unwrapping on each block. The deformation field reconstruction module is used to obtain the absolute deformation projection values ​​of the global navigation satellite system stations in each block, calculate the reference intercept of the corresponding block by combining the unwrapped phase of each block, and reconstruct the continuous absolute deformation field of the target mining area based on the reference intercept.

[0009] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. The physical property parameters of the land cover medium are retrieved by GNSS multipath reflection signals, and a medium phase delay compensation model is established based on the first principle of electromagnetic wave propagation. This effectively eliminates the systematic phase delay caused by wind and sand cover, improves the accuracy of InSAR deformation monitoring, and makes the compensation deformation results more stable and reliable.

[0010] 2. By extracting fault rupture surfaces through three-dimensional spatial clustering of underground microseismic data, generating discontinuous masks and constraining the phase unwrapping process, the cumulative propagation of unwrapping errors across faults is effectively avoided, ensuring the reliability of unwrapping results within the respective blocks on both sides of the fault. This provides a new technical approach for InSAR deformation monitoring under complex geological conditions.

[0011] 3. By determining the phase reference constant of each independent block through GNSS absolute deformation observations, the unification of multi-region references and the reconstruction of the continuous absolute deformation field of the entire mining area were realized. This fully leveraged the absolute positioning advantages of GNSS observations and the high spatial resolution advantages of InSAR observations, resulting in a high-precision absolute deformation field product with a unified spatial reference. Attached Figure Description

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

[0013] Figure 1 This is an overall block diagram of the method in Embodiment 1 of the present invention.

[0014] Figure 2 This is a flowchart illustrating the overall execution process of the method in Embodiment 1 of the present invention.

[0015] Figure 3 This is an overall block diagram of the system in Embodiment 2 of the present invention. Detailed Implementation

[0016] Hereinafter, exemplary embodiments according to this application will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are merely some embodiments of this application, and not all embodiments of this application. It should be understood that this application is not limited to the exemplary embodiments described herein.

[0017] Example 1: like Figures 1-2 As shown, a millimeter-level monitoring method for geological deformation in mining areas based on an integrated air-space-ground system includes: To facilitate understanding of the technical solution of this application, a typical application scenario is described here using an open-pit coal mine in Northwest my country. The surface of this coal mine is extensively covered with aeolian sand, ranging from approximately 0.5m to 3.0m in thickness, and frequently experiences sandstorms in winter. Furthermore, several concealed faults lie beneath the mining area. Historically, fault activation has triggered microseismic events, leading to step-deformation of the surface and severely impacting mine safety and slope stability assessments. In this embodiment, a complete monitoring period from January to December 2022 is used as an example for detailed explanation. The radar imagery utilizes C-band SAR data from the ESA Sentinel-1A satellite, with a radar wavelength of [missing information]. The incident angle is 0.056m. It is approximately 39°.

[0018] In step S1, multi-source observation data of the target mining area is acquired. This multi-source data includes time-series synthetic aperture radar (SAR) images, global navigation satellite system (GNSS) observation data, and subsurface microseismic data. The time-series SAR images refer to multiple SAR images covering the target mining area acquired within the monitoring period, used to extract surface deformation information. In this embodiment, interferometric wide-swath mode data from the European Space Agency's Sentinel-1A satellite is used, totaling 48 images spanning 12 months. The GNSS observation data includes carrier phase observations from GNSS reference stations and monitoring stations deployed in and around the mining area, used to provide high-precision absolute coordinate references and deformation references. In this embodiment, five GNSS monitoring stations are evenly deployed within the mining area, and one reference station is set up approximately 10 km outside the mining area. Each station is equipped with a dual-frequency GNSS receiver with a sampling rate of 30 seconds. Underground microseismic data refers to the catalog of mine seismic events recorded by a microseismic monitoring system installed inside the mine. It includes information such as the trigger time, magnitude, three-dimensional coordinates, and focal mechanism solution of the event. In this embodiment, 12 three-component accelerometers are arranged in multiple horizontal roadways at depths of approximately 200m to 400m underground to form a microseismic array. The monitoring system can detect microseismic events with a magnitude greater than 0.5.

[0019] It should be noted that the acquisition of the aforementioned time-series SAR images should meet coherence requirements, meaning that the time baseline of adjacent images should not be too long. Taking C-band radar as an example, it is recommended that the time baseline be controlled within 12 days, and the spatial baseline within 10% of the critical baseline, to obtain high interferometric coherence. The accuracy of GNSS observation data should be better than 5 mm horizontally and 10 mm vertically to meet the accuracy requirements of subsequent absolute deformation reconstruction. The recording of subsurface microseismic data should contain sufficient signal-to-noise ratio and event location accuracy to facilitate reliable three-dimensional cluster analysis.

[0020] In step S2, the physical property parameters of the surface covering medium of the target mining area are inverted based on the observation data of the Global Navigation Satellite System, and the phase delay of the medium is calculated based on the physical property parameters and the radar wave propagation mechanism. The original interferometric wrapping phase of the time-series synthetic aperture radar image is compensated to obtain the true geometric wrapping phase.

[0021] Traditional InSAR deformation monitoring methods typically assume radar waves propagate in free space, neglecting the influence of surface covering media on electromagnetic wave propagation speed. However, in wind-blown mining environments, sand layers, as non-magnetic media, have a significantly higher relative permittivity than air, leading to a reduction in radar wave propagation speed and thus generating additional phase delay. If this media delay is not eliminated, it will directly interfere with the deformation phase, causing systematic biases in deformation measurement results. This application addresses this problem from a physical mechanism perspective by utilizing GNSS multipath reflection signals to invert the physical properties of the medium and establishing a phase delay compensation model based on the first principles of electromagnetic wave propagation.

[0022] Specifically, the physical property parameters of the surface cover medium of the target mining area are inverted based on observation data from the Global Navigation Satellite System, including: extracting multipath reflection signals from the observation data from the Global Navigation Satellite System; and using the signal-to-noise ratio oscillation frequency and phase of the multipath reflection signals to invert the cover thickness and equivalent dielectric constant of the surface cover medium around the reference station as physical property parameters.

[0023] GNSS multipath reflection signals refer to the signals that, when a GNSS satellite signal reaches the receiver antenna, include not only the direct signal but also the reflected signal after reflection from the Earth's surface. These two signals interfere at the antenna, causing the signal-to-noise ratio (SNR) recorded by the receiver to exhibit periodic oscillations. The frequency and phase of these oscillations carry geometric and dielectric information about the Earth's surface reflecting the signal.

[0024] In this embodiment, the monitoring station numbered GNSS-03 within the mining area is selected for detailed description. This station is located on the edge of the spoil heap in the southeast of the mining area, with the surface covered by approximately 1.2m of aeolian sand. The SNR observation sequence for the GPS L1 band (frequency 1575.42MHz, wavelength approximately 0.19m) of this station from January 1st to January 31st, 2022, was extracted. By detrending the SNR time series, the residual SNR oscillation curve after removing the influence of satellite elevation angle was obtained.

[0025] Lomb-Scargle spectral analysis was performed on the residual SNR oscillation curve to extract the dominant frequency component. Based on GNSS interferometric reflectometry (GNSS-IR) theory, the SNR oscillation frequency... With antenna height Satellite elevation angle and electromagnetic wave wavelength The relationship between them is satisfied: ; By performing spectral analysis on SNR data from different elevation angle ranges, the effective height of the antenna relative to the reflecting surface, i.e., the thickness of the sand cover, can be retrieved. In this embodiment, the sand cover thickness of the GNSS-03 station is obtained through inversion. The relative error between the inversion result and the actual measured value of 1.2m is approximately 1.7%, which verifies the reliability of the inversion result.

[0026] Furthermore, the phase information of the SNR oscillation can also be used to invert the equivalent dielectric constant of the Earth's surface medium. According to Fresnel's law of reflection, when an electromagnetic wave travels at an angle of incidence... When incident on the air-sand interface, the reflection coefficient Relative permittivity of the medium Related. For horizontally polarized waves, the reflection coefficient is: ; Phase difference between reflected and direct signals It can be represented as: ; in, Let be the phase angle of the reflection coefficient. By fitting the observed SNR phase to the theoretical model, the following can be obtained: In this embodiment, the equivalent dielectric constant is obtained by inverting data from the GNSS-03 station during January 2022. This value is consistent with the typical dielectric constant range (2.5 to 3.5) for dry sand.

[0027] It should be noted that the inversion result of the equivalent dielectric constant is affected by the moisture content of the sand. When the moisture content of the sand increases, the dielectric constant increases significantly, and the attenuation of electromagnetic waves also increases. Therefore, in practical applications, meteorological data should be used to mark rainfall or snowfall events, and the quality control of the inversion results for the corresponding periods should be performed. In this embodiment, by analyzing the precipitation records of the meteorological station in the mining area, it was found that there was no effective precipitation in January 2022, and the sand was in a dry state. The inverted dielectric constant can be regarded as the case where the real part dominates, and the imaginary part (attenuation term) can be ignored.

[0028] For other GNSS stations within the mining area, the same method was used to invert the sand cover thickness and equivalent dielectric constant of each station. The inversion results showed that the sand cover thickness at the five monitoring stations varied between 0.8m and 1.5m, and the equivalent dielectric constant ranged from 2.6 to 3.2, exhibiting a certain degree of spatial heterogeneity. To extend the inversion results of these discrete stations to the entire mining area, the Kriging spatial interpolation method was used, combined with topographic elevation data and surface roughness information, to generate a sand cover thickness raster map and dielectric constant raster map covering the entire mining area, with a spatial resolution consistent with the SAR imagery (approximately 20m × 20m).

[0029] After obtaining the physical property parameters, the phase delay of the medium is calculated based on the physical property parameters and the radar wave propagation mechanism. This includes: calculating the one-way equivalent optical path difference generated by the radar wave penetrating the surface covering medium based on the principles of geometric optics and Snell's law, combined with the cover thickness, equivalent dielectric constant and radar wave incident angle; and converting the one-way equivalent optical path difference into a two-way radar interference phase delay as the medium phase delay.

[0030] According to Maxwell's equations, the speed of electromagnetic waves in a medium... Relative permittivity of the medium and relative permeability Relevant. For non-magnetic media such as sand, Therefore, the wave speed decreases to: ; in, The speed of light in a vacuum. When the radar wave travels at an angle of incidence... Passing through a thickness of When the sand layer is exposed, according to Snell's law, the angle of refraction is... satisfy: ; The actual propagation path length of radar waves in sand is The corresponding optical path length (i.e., the equivalent propagation distance in air) is Without sand cover, radar waves propagate through the air at the same vertical thickness. The corresponding slant distance path length is The difference between the two is the one-way equivalent optical path difference. : ; In this embodiment, the calculation is performed using a SAR image pixel corresponding to the GNSS-03 station as an example. The sand thickness of this pixel... m, equivalent dielectric constant Radar incident angle Substitute into the formula to calculate: ; This means that, due to the presence of sand cover, the one-way propagation path of radar waves increases by approximately 0.911m of equivalent optical path compared to the case without sand cover.

[0031] The single-path equivalent optical path difference is converted into a two-path radar interferometric phase delay as the medium phase delay. Since InSAR measures the two-way propagation phase difference of radar waves, the single-path optical path difference needs to be multiplied by 2 and converted into a phase delay. (Two-way radar interferometric phase delay) for: ; The negative sign indicates that the phase lag is caused by the medium delay. In this embodiment, the radar wavelength... m, substituted into the aforementioned calculation, yields m: ; The phase delay is attributed to Within the interval, the package phase delay is approximately rad. This delay is equivalent to approximately 6.3 mm of line-of-sight deformation error, which, if not compensated for, will directly contaminate the deformation measurement results.

[0032] For all SAR image pixels within the entire mining area, the corresponding sand thickness and dielectric constant are substituted into each pixel to calculate the dielectric phase delay of each pixel, generating a dielectric phase delay map. In this embodiment, the dielectric phase delay within the mining area is... rad to The variation between rad corresponds to an equivalent deformation error between 5 mm and 8 mm. This error level is comparable to the actual deformation signal in the mining area (usually at the centimeter level), so compensation is necessary.

[0033] The original interferometric wrapping phase of the time-series synthetic aperture radar image is compensated to obtain the true geometric wrapping phase. This includes subtracting the two-way radar interferometric phase delay from the original interferometric wrapping phase pixel by pixel to eliminate the influence of medium delay and obtain the true geometric wrapping phase.

[0034] Original interference wrap phase This is obtained by interferometric processing of temporal SAR imagery, which includes multiple components such as surface deformation phase, topographic residual phase, atmospheric delay phase, and medium delay phase. In this step, the medium phase delay is subtracted. This can eliminate the influence of dielectric delay and obtain the true geometrically enclosed phase. : ; It should be noted that, due to and All are wrap-around phases, with values ​​ranging from... Since the interval is within the specified range, after performing the phase subtraction operation, the result needs to be re-wrapped to ensure it remains within the specified range. Within the specified range. Package operations can be implemented using the following formula: ; In this embodiment, an interferometric pair consisting of two SAR images from January 1st and January 13th, 2022, is selected for illustration. Conventional interferometric processing is performed on this interferometric pair, including registration, interferogram generation, removal of flat phase, and terrain phase removal, to obtain the original interferometric wrapped phase map. At the pixel corresponding to the GNSS-03 station, the original wrapped phase... rad. Based on the aforementioned calculations, the medium phase delay of this pixel is... rad. Perform compensation calculation: ; because No need to re-wrap, true geometric wrap phase. rad.

[0035] The aforementioned compensation operation was performed on each pixel within the entire mining area to generate a compensated true geometric wrapping phase map. By comparing the wrapping phase maps before and after compensation, it can be observed that in areas with thicker sand cover (such as spoil heaps and the edges of mining pits), the phase difference before and after compensation is more significant, while in areas with exposed bedrock or thinner sand cover, the phase difference before and after compensation is smaller. This phenomenon verifies the effectiveness and necessity of medium phase delay compensation.

[0036] To further verify the compensation effect, the compensated true geometrically wrapped phase was compared with the measured deformation of the GNSS station. In this embodiment, the measured line-of-sight deformation of the GNSS-03 station from January 1st to January 13th, 2022 was... mm (the negative sign indicates the direction away from the satellite). Wrap the phase with the compensated true geometry. Convert rad to shape variable: ; This result is consistent with GNSS measured values. There is a difference of approximately 4.2 mm, which mainly stems from atmospheric delay, terrain residuals, and GNSS measurement errors. If no medium phase delay compensation is performed, the original wrapped phase can be used directly. Calculate the deformation in rad: ; The difference from the GNSS measured value is only 1.4 mm, seemingly closer. However, this closerness is accidental because the original envelope phase is mixed with the medium delay phase, causing the deformation to be underestimated. In other time phases or other regions, this accidental closeness will not be reproduced, but will instead lead to the accumulation of systematic biases. By compensating for the medium phase delay, although the difference in a single comparison may increase, from a long-term and full-area perspective, the compensated result is more stable and reliable, avoiding systematic errors caused by the variation of medium delay with time and space.

[0037] It should be noted that the accuracy of medium phase delay compensation depends on the accuracy of the physical property parameter inversion. If there are significant errors in the sand thickness or dielectric constant obtained from GNSS-IR inversion, the compensation effect will be directly affected. Therefore, in practical applications, quality control should be implemented on the inversion results, eliminating station data with excessively low signal-to-noise ratios or excessively large inversion residuals, and cross-validating them with field measurements or remote sensing imagery. Furthermore, for sand with high moisture content or significant attenuation, the imaginary part of the dielectric constant should be considered, or the corresponding temporal data should be downweighted.

[0038] The above scheme successfully enabled the inversion of physical properties of land cover media based on GNSS multipath reflection signals. Furthermore, a media phase delay compensation model was established based on the first principles of electromagnetic wave propagation, eliminating the influence of media delay from the original interferometric wrapped phase and obtaining the true geometrically wrapped phase. This design lays a high-quality data foundation for subsequent phase unwrapping and deformation reconstruction, significantly improving the accuracy and reliability of InSAR deformation monitoring in wind-blown sand-covered mining areas.

[0039] In step S3, spatial clustering is performed based on underground microseismic data to extract fault rupture surfaces and generate discontinuous masks. The discontinuous masks are superimposed on the real geometric wrapping phase to block the phase unwrapping integration path across faults, dividing the target mining area into multiple topologically isolated blocks, and performing independent phase unwrapping on each block.

[0040] Traditional InSAR phase unwrapping algorithms are based on the Nyquist sampling theorem, assuming that the phase difference between adjacent pixels is less than 1 / 2. This allows for the recovery of the true phase via an integral path. However, in mining environments with activated hidden faults, step deformation occurs on the surface on both sides of the fault, resulting in a phase difference between adjacent pixels across the fault that is much greater than that between the two faults. This severely violates the continuity assumption of phase unwrapping. Without intervention, the unwrapping algorithm will generate erroneous integer-cycle jumps at fault locations, accumulating and propagating errors across the entire map, leading to widespread systematic deviations in the deformation field. This application utilizes the three-dimensional spatial distribution of downhole microseismic events to identify the geometric morphology of the fault rupture surface from a physical mechanism perspective, projecting it onto the surface to generate a discontinuous mask. This explicitly cuts off the unwrapping integration path across faults, thereby avoiding error accumulation and ensuring the reliability of the unwrapping results within each block.

[0041] Specifically, spatial clustering based on underground microseismic data is used to extract fault rupture surfaces and generate discontinuous masks, including: performing density spatial clustering on the three-dimensional spatial coordinates contained in the underground microseismic data to extract the strike and dip angle of the main rupture surface; projecting the main rupture surface upwards onto the ground surface to generate a mask matrix with a preset buffer width as a discontinuous mask.

[0042] In this embodiment, the underground microseismic monitoring system in the mining area recorded 327 valid microseismic events from January to December 2022, with magnitudes ranging from 0.5 to 2.3. The three-dimensional coordinate distribution of these microseismic events exhibited obvious spatial clustering characteristics, reflecting the activity patterns of underground faults. To extract continuous fault rupture surfaces from the discrete microseismic event point cloud, the density-based spatial clustering algorithm DBSCAN was used for three-dimensional clustering analysis.

[0043] The core idea of ​​the DBSCAN algorithm is to identify sufficiently high-density regions in space as clusters, while treating low-density isolated points as noise points. This algorithm requires setting two key parameters: neighborhood radius. and minimum points For a specific microseismic event point If its The neighborhood contains at least If there are 1 point, then Defined as the core point; if Not a core point but located at a core point Within the neighborhood, Defined as a boundary point; otherwise These are defined as noise points. All mutually density-reachable core points and their boundary points form a cluster.

[0044] In this embodiment, a neighborhood radius is set based on geological data of the mining area and historical microseismic activity characteristics. m, minimum number of points The three-dimensional coordinates of 327 microseismic events Input the DBSCAN algorithm, where and Using planar coordinates (adopting the local coordinate system of the mining area), The coordinates are depth coordinates (zero point at the surface, positive downwards). Clustering results show that 298 of the 327 events were divided into 3 main clusters, corresponding to 3 active faults, while 29 events were identified as noise points.

[0045] Taking cluster 1 as an example, this study details 142 microseismic events, exhibiting a distinct planar spatial distribution. Principal component analysis (PCA) is used to fit the geometric parameters of the fault rupture surface from these discrete points. First, the centroid coordinates of all event points in the cluster are calculated. : ; in This represents the number of events in the cluster. The centroid coordinates are calculated as follows: m.

[0046] Then construct the covariance matrix. : ; in For the first The coordinate vector of each event Let be the centroid coordinate vector. Eigenvalue decomposition of the covariance matrix yields three eigenvalues. and its corresponding eigenvectors Minimum eigenvalue corresponding feature vector This is the normal vector of the fault rupture surface.

[0047] In this embodiment, the calculation is obtained The strike and dip angle of the fault plane can be determined based on the normal vector. Strike angle Defined as the angle between the fault strike line and true north, it can be calculated using the horizontal component of the normal vector: ; in and They are the normal vectors and Components. Substitute the numerical values ​​to obtain the orientation angle. That is, the fault strikes at 30.8° northeast.

[0048] inclination Defined as the angle between the fault plane and the horizontal plane, it can be calculated using the vertical component of the normal vector: ; in For the normal vector Components. Substitute the numerical values ​​to obtain the tilt angle. That is, the fault plane dips approximately 48.1° in the southeast direction.

[0049] Based on the above analysis, the geometric parameters of the fault corresponding to cluster 1 were successfully extracted: centroid location. The fault strikes at 30.8° NES and dips at 48.1°. The same method was used to extract the geometric parameters of the corresponding faults for clusters 2 and 3. The fault corresponding to cluster 2 strikes at 65° NES and dips at approximately 55°; the fault corresponding to cluster 3 strikes near-east-west and dips at approximately 70°.

[0050] After extracting the geometric parameters of the fault rupture surface, the main rupture surface is projected upwards onto the Earth's surface to generate a mask matrix with a preset buffer width as a discontinuous mask. Upward projection refers to extending the underground rupture surface to the Earth's surface along the dip of the fault plane to obtain the fault's trace location on the surface.

[0051] Taking the fault corresponding to cluster 1 as an example, its centroid depth is known to be 285.7 m, its dip angle is 48.1°, and its dip direction is southeast. Extending upwards along the dip direction to the surface... The projection distance is: ; Therefore, the fault's surface trace is offset approximately 257.3 m southeast relative to the horizontal projection of the centroid. Based on the strike angle of 30.8° and the dip direction, the specific coordinates of the surface trace can be calculated. In the mining area coordinate system, the fault's surface trace is a straight line striking northeast at 30.8°, passing through the coordinate point. m.

[0052] To account for the width of the fault rupture zone and the sensitivity of the phase unwrapping algorithm to boundaries, a preset buffer width is set on both sides of the surface trace to generate a mask matrix. The selection of the buffer width should comprehensively consider the actual width of the fault rupture zone, the spatial resolution of the SAR image, and the variation characteristics of the phase gradient. In this embodiment, based on geological survey data of the mining area, the width of the rupture zone of this type of fault is typically between 10m and 30m. Considering that the ground distance resolution of the Sentinel-1A image is approximately 20m, to ensure that the mask can completely cover the fault-affected area, the buffer width is set to 40m, that is, extending 40m on each side of the surface trace.

[0053] The specific steps for generating the mask matrix are as follows: First, in the geographic coordinate system of the SAR image, calculate the vertical distance from the center point of each pixel to the trace based on the direction and location of the fault surface trace. Then, determine If the pixel's width is less than the buffer width of 40m, the mask value is set to 0 (indicating a non-contiguous region); otherwise, it is set to 1 (indicating a contiguous region). This results in a binary mask matrix with the same size as the SAR image. ,in Indicates the first The cell is located within the fault discontinuity zone. This indicates that the pixel is located within a continuous region.

[0054] In this embodiment, mask matrices are generated for the three main faults, and then a logical AND operation is performed to obtain a comprehensive discontinuous mask. This mask identifies three discontinuous zones approximately 80m wide (40m on each side) within the mining area, covering a total area of ​​about 8% of the total mining area. By overlaying the mask onto the SAR image, the correspondence between the spatial distribution of fault traces and surface deformation characteristics can be visually observed.

[0055] It should be noted that the above method for extracting fault planes based on microseismic clustering assumes that the spatial distribution of microseismic events can represent the geometry of the fault. This assumption is reasonable in most cases because microseismic events are a direct manifestation of fault activity, and their spatial clustering characteristics reflect the fault rupture process. However, in some cases, microseismic events may be affected by factors such as the complexity of the stress field and the heterogeneity of the rock mass, leading to deviations in spatial distribution. Therefore, in practical applications, the clustering results should be cross-validated by combining geological data, historical seismic records, and surface deformation characteristics. If necessary, the geometric parameters or buffer width of the fault plane can be manually adjusted to ensure the accuracy of the mask.

[0056] After generating the discontinuous mask, the discontinuous mask is superimposed on the real geometry to wrap the phase, thereby blocking the phase unwrapping integration path across the fault. The target mining area is divided into multiple topologically isolated blocks, and each block is subjected to independent phase unwrapping, including: using the discontinuous mask to cut off the phase difference calculation path of adjacent pixels crossing the fault trace; and using the minimum cost flow algorithm to perform local phase unwrapping operations on each topologically isolated block segmented by the discontinuous mask.

[0057] Phase unwrapping is a crucial step in InSAR data processing, aiming to recover the true continuous phase from the wrapped phase. The range of values ​​for the wrapped phase is limited to... Within the range, the true phase may be much larger. Therefore, it is necessary to determine how many [items] should be added to each cell using an unwrapping algorithm. The phase unwrapping problem is solved by finding integer multiples of the phase cost, thus restoring the true phase. Traditional global unwrapping algorithms (such as the minimum cost flow algorithm MCF) transform the phase unwrapping problem into a network flow optimization problem by constructing a network of connections between pixels, seeking the integer periodic jump distribution that minimizes the global cost.

[0058] However, when discontinuous features such as faults exist, the true phase difference between adjacent pixels across a fault may be much greater than... This leads to the unwrapping algorithm's inability to correctly determine the location and number of cycle jumps, resulting in erroneous unwrapping results. Once this error occurs, it accumulates and propagates along the integration path, contaminating a large area of ​​the deformation field. This application introduces a discontinuity mask to explicitly cut off the integration path across faults, decomposing the global unwrapping problem into multiple local unwrapping problems, fundamentally avoiding the cross-regional propagation of errors.

[0059] Specifically, using a discontinuous mask to cut off the phase difference calculation path between adjacent cells crossing the fault trace means that when constructing the cell connection network for the phase unwrapping algorithm, for cells with a mask value of 0 (located within the discontinuous band), no connection edges are established between them and their adjacent cells, thus blocking the integration path across the fault. In the minimum cost flow algorithm, the cell connection network typically adopts a Delaunay triangulation or regular mesh structure, with each connection edge corresponding to a pair of phase differences between adjacent cells. Through the constraint of the mask, there are no longer direct connection edges between cells located on both sides of the fault. Therefore, the unwrapping algorithm cannot transmit integer-cycle jump information between the two sides of the fault, ensuring that the unwrapping processes of the two blocks are independent of each other.

[0060] In this embodiment, the true geometrically wrapped phase map obtained in step S2 is superimposed with the discontinuous mask. First, based on the mask matrix... The pixels within the mining area are divided into continuous and discontinuous regions. For mask values... For the cells that are not valid, their phase values ​​are set to invalid values ​​(such as NaN), indicating that the cell will not participate in the unwrapping calculation. Then, a Delaunay triangulation is constructed for the remaining valid cells. During the construction process, if there are cells with a mask value of 0 between the two endpoints of a certain triangle edge, the edge is deleted, thereby cutting off the connection path across the fault.

[0061] After the above processing, the originally connected mining area was divided into four topologically isolated blocks by three fault mask zones. These four blocks are spatially independent and there are no direct cell-connecting edges between them. Block 1 is located in the northern part of the mining area, accounting for approximately 35% of the total area; Block 2 is located in the central part of the mining area, accounting for approximately 30%; Block 3 is located in the southern part of the mining area, accounting for approximately 25%; and Block 4 is a small area in the southeastern corner of the mining area, accounting for approximately 10%.

[0062] After completing the block partitioning, the minimum cost flow algorithm is used to perform local phase unwrapping operations on each topologically isolated block segmented by the discontinuous mask. The basic principle of the minimum cost flow algorithm is to model the phase unwrapping problem as a network flow optimization problem: each cell is regarded as a node in the network, the connection edges between adjacent cells carry phase difference information, and the goal of unwrapping is to assign an integer-cycle jump number to each cell. This causes the unwrapping phase difference between adjacent pixels. To get as close as possible to the true phase difference while satisfying global consistency constraints.

[0063] Specifically, for adjacent pixels and Its encapsulation phase difference is The phase difference after untangling is ,in This represents the number of jumps in an integer cycle. The minimum cost flow algorithm minimizes the global cost function: ; To determine the optimal integer periodic jump distribution, where Let be the set of all connecting edges. For the edge The weights are typically set based on phase quality or coherence. This optimization problem can be solved efficiently using network flow algorithms.

[0064] In this embodiment, minimum-cost flow unwrapping is performed on each of the four blocks. Taking block 1 as an example, this block contains approximately 12,000 valid pixels, and the constructed Delaunay triangulation contains approximately 35,000 connecting edges. Based on the interferometric coherence graph, a weight is assigned to each edge: edges with high coherence have a larger weight, and edges with low coherence have a smaller weight, thereby guiding the unwrapping algorithm to propagate along high-quality paths first. After approximately 15 seconds of computation, the unwrapped phase graph of block 1 is obtained. The same method is used for blocks 2, 3, and 4 to obtain their respective unwrapped phase graphs.

[0065] It should be noted that, since the unwrapping processes of each block are independent, the unwrapped phase of each block has an arbitrary constant offset (i.e., the phase reference is not unified). This is because phase unwrapping can only determine the relative phase difference, not the absolute phase value. In traditional global unwrapping methods, a reference point is usually selected and its phase value is fixed at 0 to determine the phase reference of the entire map. However, in the partitioned unwrapping scheme of this application, the blocks are isolated by fault masks, making it impossible to unify the reference through a single reference point. Therefore, in subsequent steps, GNSS observation data is needed to determine the absolute phase reference for each block separately to achieve a unified reconstruction of the deformation field of the entire mining area.

[0066] To verify the effectiveness of the partitioned unwrapping scheme, the partitioned unwrapping results of this application are compared with those of the traditional global unwrapping method. In the traditional method, no discontinuity mask is used; the entire mining area is directly unwrapped. The comparison results show that near the fault location, the unwrapping phase of the traditional method exhibits obvious abrupt changes and stripe distortion, indicating that the unwrapping algorithm makes incorrect integer period judgments at the fault location and propagates the error to the surrounding area. In contrast, the partitioned unwrapping results of this application show a smooth and continuous phase distribution within each block, without any abnormal abrupt changes, thus verifying the effectiveness of the discontinuity mask.

[0067] Furthermore, by comparing the measured deformation with that of GNSS stations, the unwrapping accuracy can be quantitatively evaluated. In this embodiment, the five GNSS stations are located in four different blocks. For GNSS-01 station located in block 1, its measured line-of-sight deformation from January 1st to January 13th, 2022 is... mm. The unwrapped phase of the corresponding pixel at this station is converted into a deformation, resulting in... The relative error was approximately 5.9% (mm). Similar comparisons were made with four other GNSS stations, and the relative errors were all within 10%, indicating that the partitioned unwrapping results have high accuracy.

[0068] It should be noted that the successful implementation of the zonal unwrapping scheme depends on the accuracy of the discontinuous mask. If the mask position deviates from the actual fault trace, or if the buffer width is set improperly, some cross-fault connection edges may not be severed, resulting in unwrapping errors. Therefore, in practical applications, the mask parameters should be iteratively optimized by combining geological data, surface deformation characteristics, and quality assessment of the unwrapping results. If necessary, the mask boundaries can be manually edited to ensure the reliability of zonal unwrapping.

[0069] Furthermore, for situations where the number of microseismic events is small or their spatial distribution is too discrete, DBSCAN clustering may not be effective in identifying fault planes. In such cases, it is advisable to combine other geophysical exploration methods (such as seismic exploration, electrical resistivity tomography, etc.) or geological survey data to artificially pinpoint fault locations and generate masks. The core idea of ​​this application is to explicitly constrain the phase unwrapping process using discontinuous masks, and the method of mask generation can be flexibly selected according to actual data conditions, not limited to microseismic clustering methods.

[0070] The above scheme successfully achieved fault rupture surface extraction and discontinuous mask generation based on underground microseismic data. Furthermore, mask constraints enabled partitioned phase unwrapping, effectively preventing the cumulative propagation of unwrapping errors across faults and laying a reliable foundation for subsequent absolute deformation reconstruction. This design fully demonstrates the advantages of deep integration of geophysical monitoring data and remote sensing data, providing a new technical approach for deformation monitoring in mining areas under complex geological conditions.

[0071] In step S4, the absolute deformation projection values ​​of the global navigation satellite system stations in each block are obtained, the reference intercept of the corresponding block is calculated by combining the unwrapped phase of each block, and the continuous absolute deformation field of the target mining area is reconstructed based on the reference intercept.

[0072] This step is crucial for achieving unified reconstruction of the absolute deformation field across the entire mining area. In step S3, the target mining area is divided into multiple topologically isolated blocks using a discontinuous mask, and each block undergoes independent phase unwrapping. While this effectively avoids the cumulative propagation of unwrapping errors across faults, it also introduces a new problem: the unwrapped phase of each block has its own independent arbitrary constant offset, meaning the phase reference is not unified across blocks. This is because the phase unwrapping algorithm can only determine the relative phase difference, not the absolute phase value. Traditional methods typically unify the overall map reference by selecting a reference point and fixing its phase value to zero. However, in the case of partitioned unwrapping, this method cannot achieve reference transfer across blocks. This application utilizes high-precision absolute deformation observations from GNSS stations distributed within each block to determine the phase reference constant for each block, thereby achieving unified reconstruction of the deformation field across the entire mining area.

[0073] Specifically, obtaining the absolute deformation projection values ​​of the Global Navigation Satellite System (GNSS) stations within each block includes: obtaining the three-dimensional absolute deformation observation values ​​of the GNSS stations distributed within each topologically isolated block; and combining the orbital parameters and observation geometry of the radar satellites to project the three-dimensional absolute deformation observation values ​​onto the radar line of sight to obtain the absolute deformation projection values.

[0074] GNSS observation data has unique advantages: on the one hand, GNSS can provide three-dimensional absolute coordinates and deformation information, unaffected by the phase unwrapping problem; on the other hand, the location of GNSS stations is precisely known, allowing for a one-to-one correspondence with SAR image pixels. Therefore, GNSS observations can serve as an absolute reference for calibrating the constants of the InSAR unwrapped phase.

[0075] In this embodiment, a total of 5 GNSS monitoring stations are deployed within the mining area. After the partitioning and unwrapping in step S3, these 5 stations are located in 4 different blocks. The specific distribution is as follows: GNSS-01 station is located in block 1 (northern region); GNSS-02 and GNSS-03 stations are both located in block 2 (central region); GNSS-04 station is located in block 3 (southern region); and GNSS-05 station is located in block 4 (southeast corner region).

[0076] The GNSS-01 station will be used as an example for detailed explanation. From January 1st to January 13th, 2022 (corresponding to the time baseline of the first interferometric pair), the station obtained three-dimensional deformation observations through carrier phase positioning calculations. Let the components of the three-dimensional deformation vector of this station in the east, north, and vertical directions be respectively... In this embodiment, the three-dimensional deformation observation values ​​of the GNSS-01 station are: ; Among them, a positive eastward component indicates eastward movement, a positive northward component indicates northward movement, a positive vertical component indicates upward lifting, and a negative value indicates sinking.

[0077] To project GNSS 3D deformations onto the radar line-of-sight, it is necessary to combine the radar satellite's orbital parameters and observation geometry. For the Sentinel-1A satellite's interferometric wide-swath mode, its orbital parameters and observation geometry can be extracted from the satellite's precise orbital data and image metadata. Key parameters include: the satellite's 3D position vector at the time of imaging. Velocity vector and the three-dimensional coordinates of ground pixels. .

[0078] Radar line-of-sight unit vector Defined as a unit vector pointing from a ground pixel to a satellite: ; In this embodiment, for the pixel corresponding to the GNSS-01 station, the components of the radar line-of-sight unit vector in the east, north, and vertical directions are calculated based on satellite orbit data as follows: ; These components reflect the satellite's observation geometry: a positive eastward component indicates the satellite is east of the target pixel, while a negative vertical component indicates the satellite is above the target pixel (for ascending orbits). Line of sight to angle of incidence. It can be calculated using the vertical component: ; This angle matches the incident angle of 39° given in step S1, verifying the correctness of the observation geometry parameters.

[0079] The GNSS three-dimensional deformation observations are projected onto the radar line-of-sight to obtain the absolute deformation projection values. : ; Substitute the values ​​from the GNSS-01 station: ; It should be noted that the sign convention used in the above projection calculations is: positive A negative value indicates that the ground target is moving away from the satellite along the line of sight (typically corresponding to ground subsidence or northward displacement), while a negative value indicates movement towards the satellite. In this embodiment, the absolute deformation projection value of the GNSS-01 station is... mm indicates that the station moved approximately 8.9 mm away from the satellite along the line of sight during the observation period, which is consistent with the land subsidence trend in the region.

[0080] The same method was used to perform projection calculations on other GNSS stations within the mining area to obtain the absolute deformation projection values ​​for each station: ; ; ; ; The data above shows that GNSS stations in different blocks exhibit differentiated deformation characteristics: Blocks 1 and 2 show deformation away from the satellite (positive values), Block 3 shows deformation towards the satellite (negative values), and Block 4 shows slight deformation away from the satellite. This differentiated deformation pattern is a typical characteristic of fault activity: the surface deformation on both sides of the fault is in opposite directions or the magnitudes differ significantly.

[0081] The reference intercept of each block is calculated by combining the unwrapped phase of each block, including: extracting the unwrapped phase of the block corresponding to the pixel location of the Global Navigation Satellite System station; and calculating the phase reference constant of each block as the reference intercept based on the difference between the phase equivalent corresponding to the absolute deformation projection value and the unwrapped phase of the block.

[0082] The physical meaning of the phase reference constant is the offset between the unwrapped phase and the absolute phase. For the first... GNSS stations within each block Its unwrapping phase in SAR imagery is The absolute deformation projection value measured by GNSS The corresponding phase equivalent should be The difference between the two is the phase reference constant for that block. : ; The derivation of this formula is based on the following physical principle: there is a definite relationship between the InSAR interferometric phase and the line-of-sight deformation. For two-way radar interferometry, the deformation... The resulting phase change is (The negative sign indicates that the deformation farther from the satellite corresponds to a phase lag). Therefore, if the absolute deformation measured by GNSS is known... Then the absolute phase that the station should have can be calculated: ; InSAR unwrapping phase with absolute phase There is a constant offset between them: ; Therefore, the phase reference constant is: ; It should be noted that the symbol conventions in the above formulas should be consistent with the actual data processing. In this embodiment, the GNSS deformation projection value... A positive value indicates distance from the satellite, corresponding to a negative phase lag; therefore, the absolute phase is negative. InSAR unwrapped phase is typically expressed in radians, with a positive value indicating a phase lead. In specific calculations, a consistent notation should be used to avoid confusion. For simplicity, this embodiment uses the formula form from the technical disclosure, defining the phase reference constant as: ; Under this definition, a positive reference constant indicates that the unwrapped phase has a negative offset relative to the GNSS absolute reference, and this constant needs to be added for correction.

[0083] In this embodiment, the calculation is performed using block 1 as an example. The GNSS-01 station is located in block 1, and its absolute deformation projection value... Extract the unwrapped phase value of the corresponding pixel of the station from the unwrapped phase map of block 1 obtained in step S3. Radar wavelength Substitute into the formula to calculate: ; Therefore, the phase reference constant of block 1 This result indicates that the unwrapping phase of block 1 has a positive offset of approximately 0.153 rad relative to the absolute reference, which needs to be subtracted (or equivalently added) in subsequent reconstruction. (The reference constant in rad).

[0084] The same method is used to calculate the phase reference constants for other blocks. For block 2, since there are two GNSS stations (GNSS-02 and GNSS-03), the two reference constants can be calculated separately and then averaged to improve robustness. ; ; ; For Block 3 and Block 4, there is only one GNSS station each: ; ; Through the above calculations, the phase reference constants for each block were obtained: , , , These constants will be used for subsequent reconstruction of the absolute deformation field.

[0085] It should be noted that the absolute differences in the reference constants reflect the relative shifts in the unwrapping phases of each block. The differences in reference constants between different blocks can be significant, which directly reflects the inconsistency in references caused by partitioned unwrapping. By determining the reference constants for each block, a unified reconstruction of the deformation field across the entire mining area can be achieved.

[0086] The continuous absolute deformation field of the target mining area is reconstructed based on the reference intercept, including: summing the unwrapped phase of any pixel in each block with the phase reference constant of the corresponding block; converting the summed absolute phase into the physical deformation of the radar line of sight direction, and splicing them together to reconstruct the continuous absolute deformation field of the unified spatial reference of the entire mining area.

[0087] For the Any pixel within a block Its absolute deformation The calculation formula is: ; The physical meaning of this formula is: adding a reference constant to the unwrapped phase yields the absolute phase with reference to the GNSS absolute coordinate system; then, the absolute phase is converted into physical deformations. (Coefficients) It is a scaling factor that converts phase (radians) into deformation (meters).

[0088] In this embodiment, a typical cell within block 1 is used as an example for calculation. Let the unwrapping phase of this cell be... The baseline constant of block 1 Substitute into the formula: ; The absolute deformation of the pixel mm indicates that the satellite moved approximately 15.0 mm away from the satellite along the line of sight during the observation period.

[0089] The above calculations are performed on all valid pixels within block 1 to obtain the absolute deformation field of block 1. Similarly, calculations are performed on blocks 2, 3, and 4 respectively. Taking block 2 as an example, let the unwrapping phase of a certain pixel be... reference constant : ; For block 3, let the untangling phase of a certain cell be... reference constant : ; For block 4, let the untangling phase of a certain cell be... reference constant : ; After calculating the absolute deformation of each block, the deformation fields of each block are stitched together to reconstruct a continuous absolute deformation field with a unified spatial reference for the entire mining area. During the stitching process, the deformation values ​​at the boundaries of each block should remain consistent. However, since the fault mask area is excluded from the calculation, deformation discontinuities will appear at the fault traces. This discontinuity reflects a real physical phenomenon, rather than a calculation error.

[0090] In this embodiment, the reconstructed absolute deformation field of the mining area exhibits the following characteristics: the deformation of Block 1 (northern region) is between +5mm and +25mm, showing an overall subsidence trend; the deformation of Block 2 (central region) is between +2mm and +18mm, with a subsidence amplitude slightly smaller than that of Block 1; the deformation of Block 3 (southern region) is between -8mm and +12mm, with some areas showing an uplift trend; and the deformation of Block 4 (southeast corner region) is between +3mm and +15mm, showing slight subsidence.

[0091] Spatially, the deformation fields on both sides of the fault exhibit a clear step-like characteristic. Taking the fault dividing Block 1 and Block 2 as an example, the fault strikes 30.8°E. The deformation on the northwest side (Block 1) is generally greater than that on the southeast side (Block 2), with a difference in deformation between 5 mm and 10 mm. This step-like deformation pattern is consistent with the differential subsidence caused by fault activity, verifying the rationality and effectiveness of the zonal unwrapping scheme.

[0092] To further verify the accuracy of the absolute deformation field reconstruction, the reconstruction results were compared with the measured deformation from GNSS stations. Since the GNSS station observations have already been used to determine the reference constants, this comparison can only verify the degree of agreement between the stations, not the accuracy independently. Taking the GNSS-01 station as an example, its measured absolute deformation projection value is... The absolute deformation of the reconstructed cell is: ; The two results are completely consistent, verifying the correctness of the baseline constant calculation and absolute deformation reconstruction.

[0093] It should be noted that the multi-zone benchmark unification scheme of this application relies on the existence of at least one GNSS station within each block. If there is no GNSS station in a block, the benchmark constants for that block cannot be determined, and other methods need to be used for benchmark transfer. For example, relative benchmark alignment can be performed using topographic feature points on both sides of a fault, or a regional network adjustment method can be used to perform joint calculations using common feature points across blocks. In practical applications, the GNSS station deployment scheme should be rationally planned according to the geological conditions and monitoring needs of the mining area to ensure that there is at least one GNSS station in each potential independent block.

[0094] Furthermore, the deformation observation accuracy of GNSS stations directly affects the reconstruction accuracy of the absolute deformation field. In this embodiment, the accuracy of GNSS carrier phase positioning is better than 5mm (horizontal) and 10mm (vertical), and the accuracy projected onto the line-of-sight back is approximately 6mm to 8mm. Therefore, the theoretical accuracy of the reconstructed absolute deformation field is approximately one phase period (approximately 28mm for C-band) plus the GNSS projection error. Further improvements in the accuracy and reliability of the deformation field can be achieved through long-term time-series observations and joint calculations using multiple interferometric pairs.

[0095] The above scheme successfully achieved multi-region benchmark unification and reconstruction of the continuous absolute deformation field of the entire mining area based on GNSS absolute deformation constraints. This design fully utilizes the absolute positioning advantages of GNSS observations and the high spatial resolution advantages of InSAR observations, effectively solving the benchmark inconsistency problem caused by regional unwrapping, and providing reliable absolute deformation field products for deformation monitoring in mining areas under complex geological conditions.

[0096] In summary, the deformation monitoring method for mining areas based on aeolian dielectric compensation and fault topological constraints proposed in this application achieves the following technical effects: First, the physical property parameters of the land cover medium are retrieved by inverting GNSS multipath reflection signals, and a medium phase delay compensation model is established based on the first principle of electromagnetic wave propagation. This effectively eliminates the systematic phase delay caused by wind and sand cover, improves the accuracy of InSAR deformation monitoring, and makes the compensation deformation results more stable and reliable.

[0097] Second, by extracting fault rupture surfaces through three-dimensional spatial clustering of underground microseismic data, generating discontinuous masks and constraining the phase unwrapping process, the cumulative propagation of unwrapping errors across faults is effectively avoided, ensuring the reliability of unwrapping results within the respective blocks on both sides of the fault, and providing a new technical approach for InSAR deformation monitoring under complex geological conditions.

[0098] Third, by using GNSS absolute deformation observations to determine the phase reference constants of each independent block, the unification of multi-region references and the reconstruction of the continuous absolute deformation field of the entire mining area were achieved. This fully leveraged the absolute positioning advantages of GNSS observations and the high spatial resolution advantages of InSAR observations, resulting in high-precision absolute deformation field products with a unified spatial reference.

[0099] Fourth, this application deeply integrates multi-source remote sensing data (time-series SAR images, GNSS observation data, and underground microseismic data), forming a complete technical chain from data acquisition, physical parameter inversion, phase compensation, partition unwrapping to absolute deformation reconstruction. It provides a systematic solution for deformation monitoring in complex mining areas with wind and sand cover and fault development, and has significant practical value and prospects for promotion.

[0100] Example 2: like Figure 3 As shown, the millimeter-level monitoring system for geological deformation in mining areas, based on an integrated air-space-ground system, includes: The data acquisition module is used to acquire multi-source observation data of the target mining area, including time-series synthetic aperture radar imagery, global navigation satellite system observation data, and underground microseismic data. The phase compensation module is used to invert the physical property parameters of the surface covering medium of the target mining area based on observation data from the Global Navigation Satellite System, and calculate the phase delay of the medium based on the physical property parameters and the radar wave propagation mechanism. It compensates the original interferometric wrapping phase of the time-series synthetic aperture radar image to obtain the true geometric wrapping phase. The phase unwrapping module is used to perform spatial clustering based on underground microseismic data to extract fault rupture surfaces and generate discontinuous masks. The discontinuous masks are superimposed on the real geometric wrapping phase to block the phase unwrapping integration path across faults, dividing the target mining area into multiple topologically isolated blocks, and performing independent phase unwrapping on each block. The deformation field reconstruction module is used to obtain the absolute deformation projection values ​​of the global navigation satellite system stations in each block, calculate the reference intercept of the corresponding block by combining the unwrapped phase of each block, and reconstruct the continuous absolute deformation field of the target mining area based on the reference intercept.

[0101] The above description is merely an example and illustration of the structure of the present invention. Those skilled in the art can make various modifications or additions to the specific embodiments described, or use similar methods to replace them, as long as they do not deviate from the structure of the invention or exceed the scope defined in the claims, all of which should fall within the protection scope of the present invention.

[0102] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0103] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to any specific implementation. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A millimeter-level monitoring method for geological deformation in mining areas based on an integrated air-space-ground system, characterized in that, include: Acquire multi-source observation data of the target mining area, including time-series synthetic aperture radar imagery, global navigation satellite system observation data, and underground microseismic data; Based on the observation data of the Global Navigation Satellite System, the physical property parameters of the surface covering medium of the target mining area are inverted, and the phase delay of the medium is calculated based on the physical property parameters and the radar wave propagation mechanism. The original interferometric wrapping phase of the time-series synthetic aperture radar image is compensated to obtain the true geometric wrapping phase. Spatial clustering is performed based on the underground microseismic data to extract fault rupture surfaces and generate discontinuous masks. The discontinuous masks are superimposed on the real geometrically wrapped phase to block the phase unwrapping integration path across faults. The target mining area is divided into multiple topologically isolated blocks, and each block is subjected to independent phase unwrapping. The absolute deformation projection values ​​of the global navigation satellite system stations in each block are obtained. The reference intercept of the corresponding block is calculated by combining the unwrapped phase of each block. The continuous absolute deformation field of the target mining area is reconstructed based on the reference intercept.

2. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 1, characterized in that, The physical property parameters of the surface cover medium of the target mining area are retrieved based on the observation data of the Global Navigation Satellite System, including: Extract the multipath reflection signal from the observation data of the Global Navigation Satellite System; The signal-to-noise ratio oscillation frequency and phase of the multipath reflection signal are used to invert the coverage thickness and equivalent dielectric constant of the surface covering medium around the reference station as the physical property parameters.

3. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 2, characterized in that, The calculation of medium phase delay based on the aforementioned physical property parameters and radar wave propagation mechanism includes: Based on the principles of geometric optics and Snell's law, combined with the coverage thickness, the equivalent dielectric constant and the radar wave incident angle, the one-way equivalent optical path difference generated by radar waves penetrating the surface coverage medium is calculated. The single-path equivalent optical path difference is converted into a two-path radar interferometric phase delay as the medium phase delay.

4. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 3, characterized in that, The original interferometric wrapping phase of the time-series synthetic aperture radar image is compensated to obtain the true geometric wrapping phase, including: The true geometrically wrapped phase is obtained by subtracting the two-way radar interferometric phase delay pixel by pixel from the original interferometric wrapped phase to eliminate the influence of medium delay.

5. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 1, characterized in that, Spatial clustering based on the aforementioned subsurface microseismic data to extract fault rupture surfaces and generate discontinuity masks includes: Density space clustering is performed on the three-dimensional spatial coordinates contained in the underground microseismic data to extract the strike and dip angle of the main rupture surface; The main fracture surface is projected onto the ground surface to generate a mask matrix with a preset buffer width as the discontinuous mask.

6. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 5, characterized in that, The discontinuous mask is superimposed on the true geometrically wrapped phase to block the phase unwrapping integration path across faults, dividing the target mining area into multiple topologically isolated blocks, and performing independent phase unwrapping on each block, including: The discontinuous mask is used to cut off the phase difference calculation path of adjacent pixels crossing the fault trace; The minimum cost flow algorithm is used to perform local phase unwrapping operations on each topologically isolated block segmented by the discontinuous mask.

7. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 1, characterized in that, Obtain the absolute deformation projection values ​​of Global Navigation Satellite System (GNSS) stations within each block, including: Obtain three-dimensional absolute deformation observations of Global Navigation Satellite System (GNSS) sites distributed within each topologically isolated block; By combining the orbital parameters and observation geometry of the radar satellite, the observed three-dimensional absolute deformation values ​​are projected onto the radar line of sight to obtain the projected absolute deformation values.

8. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 7, characterized in that, The reference intercept of the corresponding block is calculated by combining the unwrapping phase of each block, including: Extract the unwrapped phase of the block corresponding to the pixel location of the Global Navigation Satellite System station; Based on the difference between the phase equivalent corresponding to the absolute deformation projection value and the unwrapping phase of the block, the phase reference constant of each block is calculated as the reference intercept.

9. The millimeter-level monitoring method for geological deformation in mining areas based on integrated air-space-ground systems as described in claim 8, characterized in that, Reconstructing the continuous absolute deformation field of the target mining area based on the aforementioned benchmark intercept includes: The unwrapped phase of any pixel within each block is summed with the phase reference constant of the corresponding block; The summed absolute phase is converted into physical deformation in the radar line-of-sight direction, and then spliced ​​together to reconstruct a continuous absolute deformation field with a unified spatial reference for the entire mining area.

10. A millimeter-level monitoring system for geological deformation in mining areas based on an integrated air-space-ground system, characterized by: The method for monitoring millimeter-level geological deformation in mining areas based on an integrated air-space-ground system as described in any one of claims 1-9 includes: The data acquisition module is used to acquire multi-source observation data of the target mining area, including time-series synthetic aperture radar imagery, global navigation satellite system observation data, and underground microseismic data. The phase compensation module is used to invert the physical property parameters of the surface covering medium of the target mining area based on the observation data of the global navigation satellite system, and calculate the phase delay of the medium based on the physical property parameters and the radar wave propagation mechanism, and compensate the original interferometric wrapping phase of the time-series synthetic aperture radar image to obtain the true geometric wrapping phase. The phase unwrapping module is used to perform spatial clustering based on the underground microseismic data to extract fault rupture surfaces and generate discontinuous masks. The discontinuous masks are superimposed on the real geometric wrapping phase to block the phase unwrapping integration path across faults, dividing the target mining area into multiple topologically isolated blocks, and performing independent phase unwrapping on each block. The deformation field reconstruction module is used to obtain the absolute deformation projection values ​​of the global navigation satellite system stations in each block, calculate the reference intercept of the corresponding block by combining the unwrapped phase of each block, and reconstruct the continuous absolute deformation field of the target mining area based on the reference intercept.

Images