Earth surface deformation monitoring method and system based on time sequence InSAR

Abstract

The invention is suitable for the technical field of earth surface monitoring, and provides an earth surface deformation monitoring method and system based on a time sequence InSAR, and the method comprises the following steps: carrying out the SAR image screening and interferogram generation, and collecting corresponding meteorological data and thermal infrared data; performing space-time adaptive atmospheric phase correction, taking the air pressure vertical gradient and the temperature anomaly as driving factors of atmospheric delay, and separating an atmospheric phase through a space-time weighted model; dynamic deformation modeling is carried out, and deformation is decomposed into linear deformation and nonlinear deformation; the method comprises the following steps: taking a mining area road network as a geometric constraint, unwrapping a coherent region by adopting a minimum cost flow algorithm, converting an unwrapping phase into sight-line-direction deformation, and calculating horizontal and vertical deformation components in combination with InSAR sight-line-direction deformation and a digital elevation model. The method adapts to a complex deformation mechanism of a mining area by capturing linear and nonlinear deformation. Through sight line deformation and DEM geometric projection, vertical and horizontal deformation separation is realized, and the deformation direction is determined.

Description

Technical Field

[0001] This invention relates to the field of surface monitoring technology, specifically to a method and system for monitoring surface deformation based on time-series InSAR. Background Technology

[0002] InSAR extracts surface deformation through phase difference in radar images, but it is limited by factors such as atmospheric delay, decoherence, and ambiguity of phase unwrapping, making it difficult to directly obtain high-precision three-dimensional deformation. Surface deformation in mining areas is complex (long-term subsidence, instantaneous blasting, slope slippage, etc.), requiring the integration of multi-source data (meteorological, thermal infrared, DEM) and geometric constraints to improve monitoring accuracy. Atmospheric phase correction relies on external data (such as the ECMWF meteorological model), but its resolution is insufficient; traditional deformation models (such as linear models) cannot capture nonlinear deformation; and phase unwrapping lacks geometric constraints, easily leading to the propagation of local errors. Therefore, a surface deformation monitoring method and system based on time-series InSAR is needed to address these problems. Summary of the Invention

[0003] In view of the shortcomings of the existing technology, the purpose of this invention is to provide a method and system for monitoring surface deformation based on time-series InSAR, so as to solve the problems existing in the above-mentioned background technology.

[0004] This invention is implemented as follows: a method for monitoring surface deformation based on time-series InSAR, the method comprising the following steps:

[0005] SAR image screening and interferogram generation are performed, and corresponding meteorological and thermal infrared data are collected. The surface temperature is retrieved through thermal infrared data.

[0006] Spatiotemporal adaptive atmospheric phase correction is performed, taking the vertical pressure gradient and temperature anomaly as the driving factors of atmospheric delay, and separating the atmospheric phase through a spatiotemporal weighted model;

[0007] Dynamic deformation modeling is performed, decomposing deformation into linear and nonlinear deformation, and capturing the deformation mechanism of the mining area through segmented modeling;

[0008] The mining area road network is used as a geometric constraint. The minimum cost flow algorithm is used to unwrap the coherent region and convert the unwrapped phase into line-of-sight deformation. The horizontal and vertical deformation components are calculated by combining InSAR line-of-sight deformation with digital elevation model.

[0009] As a further aspect of the present invention: the steps of SAR image screening and interferogram generation specifically include:

[0010] SAR images with matching spatial and temporal baselines were selected, and SAR images with cloud cover exceeding 10% were removed.

[0011] A small baseline set strategy is adopted to construct a temporal interferometric network and perform multi-view processing on the interferogram;

[0012] Calculate the coherence coefficient threshold and remove low-coherence regions.

[0013] As a further aspect of the present invention: the step of performing spatiotemporal adaptive atmospheric phase correction specifically includes:

[0014] Retrieve ERA5 data, extract PT (tropopause pressure) and PS (surface pressure) from ERA5 data, and calculate the vertical pressure gradient ΔP / Δh = (PS-PT) / H, where H is the atmospheric elevation.

[0015] Calculate the residual between surface temperature and ERA5 temperature, identify temperature anomaly areas, and determine the temperature anomaly value ΔT;

[0016] Construct a spatiotemporal weighted model and determine the spatial weight W1 = exp(-d 2 / (2σs 2 )), d is the pixel spacing, σs = 500m, and the time weight is determined as W2 = exp(-Δt) 2 / (2σt 2 )), Δt is the time difference, σt = 30 days; determine the atmospheric phase delay Φatm = Σ[WsWt(aΔP / Δh+bΔT)], and solve for the coefficients a and b using the least squares method;

[0017] The atmospheric phase is corrected, phase filtering is performed, and the corrected phase is fitted with a quadratic polynomial to eliminate residual trend terms.

[0018] As a further aspect of the present invention: the step of performing dynamic deformation modeling specifically includes:

[0019] Extract the time series phase Δφ(t,r) from the atmospherically corrected interferogram set;

[0020] To perform linear deformation parameter inversion, for each pixel r, establish a system of linear equations: Δφl(ti,r)=4π / λ×v(r)×ti+∈i, where λ is the radar wavelength, v(r) is the linear deformation rate, and ∈i represents the residual term;

[0021] Subtracting the linear deformation phase from the total phase yields the nonlinear deformation phase Δφn(t,r) = Δφ(t,r) - Δφl(t,r);

[0022] Perform piecewise nonlinear inversion, and establish a system of nonlinear equations for each pixel r: Δφn(ti,r)=4π / λ×A(r)×(1-e -ti / τ(r) )+∈i, A(r) represents the amplitude A(r), and τ(r) represents the decay time constant;

[0023] Generate deformation rate cloud maps and mark high-risk settlement areas; generate deformation amplitude cloud maps and decay time constant cloud maps.

[0024] As a further aspect of the present invention: the step of calculating the horizontal and vertical deformation components specifically includes:

[0025] Extract satellite orbital parameters from the SAR image header file. These parameters include the ascending / descending orbit mode, incident angle, and azimuth angle.

[0026] Based on the digital elevation model, the slope and aspect of each pixel are calculated. The slope represents the degree of inclination of the ground surface, and the aspect represents the direction of inclination of the ground surface.

[0027] Line-of-sight deformation projection decomposition is performed. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle, while the vertical deformation component is determined by the incident angle.

[0028] Another object of the present invention is to provide a surface deformation monitoring system based on time-series InSAR, the system comprising:

[0029] The image filtering module is used to filter SAR images and generate interferograms, collect corresponding meteorological data and thermal infrared data, and retrieve the surface temperature through thermal infrared data.

[0030] The atmospheric phase correction module is used for spatiotemporal adaptive atmospheric phase correction. It takes the vertical pressure gradient and temperature anomaly as the driving factors of atmospheric delay and separates the atmospheric phase through a spatiotemporal weighted model.

[0031] The deformation mechanism determination module is used for dynamic deformation modeling, decomposing deformation into linear and nonlinear deformation, and capturing the deformation mechanism of the mining area through segmented modeling.

[0032] The deformation component determination module is used to retrieve the road network in the mining area as a geometric constraint. It adopts the minimum cost flow algorithm to unwrap the coherent region, convert the unwrapped phase into line-of-sight deformation, and calculate the horizontal and vertical deformation components by combining the InSAR line-of-sight deformation with the digital elevation model.

[0033] As a further aspect of the present invention: the image screening module includes:

[0034] The SAR image filtering unit is used to filter out SAR images whose spatial and temporal baselines match, and to remove SAR images with cloud cover greater than 10%.

[0035] Interferometric network building unit, used to construct temporal interferometric networks using a small baseline set strategy, and to perform multi-view processing on interferograms;

[0036] The coherence coefficient calculation unit is used to calculate the coherence coefficient threshold and remove low-coherence regions.

[0037] As a further aspect of the present invention: the atmospheric phase correction module includes:

[0038] The vertical pressure gradient unit is used to retrieve ERA5 data, extract the tropopause pressure PT and the surface pressure PS from the ERA5 data, and calculate the vertical pressure gradient ΔP / Δh=(PS-PT) / H, where H is the atmospheric elevation.

[0039] The temperature anomaly unit is used to calculate the residual between the surface temperature and the ERA5 temperature, identify temperature anomaly areas, and determine the temperature anomaly value ΔT.

[0040] Spatiotemporal weighted model unit, used to construct spatiotemporal weighted model and determine spatial weight W1 = exp(-d 2 / (2σs 2 )), d is the pixel spacing, σs = 500m, and the time weight is determined as W2 = exp(-Δt) 2 / (2σt 2 )), Δt is the time difference, σt = 30 days; determine the atmospheric phase delay Φatm = Σ[WsWt(aΔP / Δh+bΔT)], and solve for the coefficients a and b using the least squares method;

[0041] The atmospheric phase correction unit is used to correct the atmospheric phase, perform phase filtering, and perform quadratic polynomial fitting on the corrected phase to eliminate residual trend terms.

[0042] As a further aspect of the present invention: the deformation mechanism determination module includes:

[0043] The time series phase extraction unit is used to extract the time series phase Δφ(t,r) from the atmospherically corrected interferogram set.

[0044] The linear deformation inversion unit is used to perform linear deformation parameter inversion. For each pixel r, a system of linear equations is established: Δφl(ti,r)=4π / λ×v(r)×ti+∈i, where λ is the radar wavelength, v(r) is the linear deformation rate, and ∈i represents the residual term;

[0045] Nonlinear deformation phase unit, used to subtract linear deformation phase from total phase to obtain nonlinear deformation phase Δφn(t,r)=Δφ(t,r)-Δφl(t,r);

[0046] The nonlinear inversion unit is used to perform piecewise nonlinear inversion. For each pixel r, a system of nonlinear equations is established: Δφn(ti,r)=4π / λ×A(r)×(1-e -ti / τ(r) )+∈i, A(r) represents the amplitude A(r), and τ(r) represents the decay time constant;

[0047] The cloud map generation unit is used to generate deformation rate cloud maps and mark high-risk settlement areas; it also generates deformation amplitude cloud maps and decay time constant cloud maps.

[0048] As a further aspect of the present invention: the deformation component determination module includes:

[0049] The orbit parameter extraction unit is used to extract satellite orbit parameters from the SAR image header file. The satellite orbit parameters include the ascending / descending orbit mode, incident angle, and azimuth angle.

[0050] The slope and aspect calculation unit is used to calculate the slope and aspect of each pixel based on the digital elevation model. The slope represents the degree of inclination of the ground surface, and the aspect represents the direction of inclination of the ground surface.

[0051] The horizontal and vertical deformation elements are used to perform line-of-sight deformation projection decomposition. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle, while the vertical deformation component is determined by the incident angle.

[0052] Compared with the prior art, the beneficial effects of the present invention are:

[0053] This invention separates atmospheric phase by synchronizing meteorological and thermal infrared data and combining them with a spatiotemporal weighted model, reducing reliance on external models. It introduces dynamic deformation modeling to capture linear and nonlinear deformations in segments, adapting to the complex deformation mechanisms of mining areas. Utilizing the mining area road network as a geometric constraint, combined with a minimum cost flow algorithm, it improves the accuracy of unwrapped phase. Through the geometric projection of line-of-sight deformation onto the DEM, it achieves separation of vertical and horizontal deformations, clearly defining the deformation direction. Attached Figure Description

[0054] Figure 1 This is a flowchart of a method for monitoring surface deformation based on time-series InSAR.

[0055] Figure 2 This is a flowchart of SAR image selection in a time-series InSAR-based surface deformation monitoring method.

[0056] Figure 3 This is a flowchart of adaptive atmospheric phase correction in a time-series InSAR-based method for monitoring surface deformation.

[0057] Figure 4 This is a flowchart of dynamic deformation modeling in a time-series InSAR-based surface deformation monitoring method.

[0058] Figure 5 This is a flowchart illustrating the calculation of deformation components in a time-series InSAR-based surface deformation monitoring method.

[0059] Figure 6This is a schematic diagram of a surface deformation monitoring system based on time-series InSAR. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0061] The specific implementation of the present invention will be described in detail below with reference to specific embodiments.

[0062] like Figure 1 As shown, this embodiment of the invention provides a method for monitoring surface deformation based on time-series InSAR, the method comprising the following steps:

[0063] S100 performs SAR image screening and interferogram generation, collects corresponding meteorological data and thermal infrared data, and retrieves surface temperature through thermal infrared data;

[0064] S200 performs spatiotemporal adaptive atmospheric phase correction, taking the vertical pressure gradient and temperature anomaly as the driving factors of atmospheric delay, and separating the atmospheric phase through a spatiotemporal weighted model;

[0065] S300 performs dynamic deformation modeling, decomposing deformation into linear and nonlinear deformation, and capturing the deformation mechanism of the mining area through segmented modeling.

[0066] S400 uses the mining area road network as a geometric constraint, employs a minimum cost flow algorithm to unwrap the coherent region, converts the unwrapped phase into line-of-sight deformation, and combines InSAR line-of-sight deformation with a digital elevation model to calculate the horizontal and vertical deformation components.

[0067] In this embodiment of the invention, SAR image pairs that meet the spatiotemporal baseline threshold are selected, and interferograms are generated to extract the surface deformation phase. Meteorological data (such as air pressure and humidity) and thermal infrared data are simultaneously acquired. The surface temperature (LST) is retrieved via thermal infrared retrieval for subsequent atmospheric phase correction. This ensures the quality of the interferograms (low baseline reduces decoherence) while simultaneously acquiring environmental parameters (temperature and air pressure) to quantify the interference of atmospheric delay on deformation measurements. Then, the vertical pressure gradient and temperature anomalies are used as the main driving factors of atmospheric delay to construct a spatiotemporal weighted model (such as Gaussian weighted or Kriging interpolation) to separate the atmospheric phase from the surface deformation phase. This eliminates the interference of atmospheric turbulence and stratification delay on deformation measurements, improving deformation inversion accuracy (error reduced to the millimeter level). Next, the deformation is decomposed into linear deformation (such as long-term subsidence) and nonlinear deformation (such as instantaneous deformation during blasting). Segmented modeling (linear least squares + nonlinear optimization) captures the complex deformation mechanism of the mining area, separating different deformation modes to avoid insufficient fitting of complex deformation by a single model and improve the analytical capability of the deformation time series. Finally, the mining area road network is used as a geometric constraint (road deformation is usually continuous and stable). The minimum cost flow algorithm (MCF) is used to unwrap the coherent region, converting the unwrapped phase into line-of-sight (LOS) deformation. Combined with a digital elevation model (DEM), the LOS deformation is decomposed into vertical, east-west, and north-south deformation components through geometric projection relationships, effectively resolving phase unwrapping ambiguities, achieving high-precision three-dimensional deformation inversion, and clarifying the deformation direction and spatial distribution.

[0068] like Figure 2 As shown, in a preferred embodiment of the present invention, the steps of SAR image screening and interferogram generation specifically include:

[0069] S101, filter out SAR images that match the spatial baseline and temporal baseline, and remove SAR images with cloud cover greater than 10%.

[0070] S102 employs a small baseline set strategy to construct a temporal interferometric network and perform multi-view processing on the interferogram;

[0071] S103, calculate the coherence coefficient threshold and remove low coherence regions.

[0072] In this embodiment of the invention, when selecting suitable SAR images, image pairs with a spatial baseline <200m and a temporal baseline <45 days are prioritized, while SAR images with cloud cover higher than 10% are removed, and clear weather data are given priority. Then, a Small Baseline Set (SBAS) strategy is employed to construct a temporal interferometric network for interferometric processing: precise registration (sub-pixel level, residual ≤0.1 pixels), multi-look processing (azimuth × range = 4 × 2), and adaptive filtering (Goldstein filtering). Finally, a coherence coefficient threshold (>0.3) is calculated, and low-coherence areas, such as vegetated areas, are removed.

[0073] like Figure 3 As shown, in a preferred embodiment of the present invention, the step of performing spatiotemporal adaptive atmospheric phase correction specifically includes:

[0074] S201, retrieve ERA5 data, extract the tropopause pressure PT and surface pressure PS from the ERA5 data, and calculate the vertical pressure gradient ΔP / Δh=(PS-PT) / H;

[0075] S202, calculate the residual between the surface temperature and the ERA5 temperature to determine the temperature anomaly value ΔT;

[0076] S203, Construct a spatiotemporal weighted model and determine the spatial weight W1 = exp(-d 2 / (2σs 2 Determine the time weight W2 = exp(-Δt) 2 / (2σt 2 Determine the atmospheric phase delay Φatm=Σ[WsWt(aΔP / Δh+bΔT)], and solve for the coefficients a and b using the least squares method;

[0077] S204, corrects the atmospheric phase, performs phase filtering, and then performs quadratic polynomial fitting on the corrected phase.

[0078] In this embodiment of the invention, when performing atmospheric phase correction, ERA5 data needs to be retrieved. The tropopause pressure PT and surface pressure PS are extracted from the ERA5 data, and the vertical pressure gradient ΔP / Δh = (PS-PT) / H is calculated, where H is the atmospheric elevation, taken as 8 kilometers. Then, the residual between the surface temperature and the ERA5 temperature is calculated to identify temperature anomaly regions and determine the temperature anomaly value ΔT. Next, a spatiotemporal weighted model is constructed, and the spatial weight W1 = exp(-d 2 / (2σs 2 )), d is the pixel spacing, σs = 500m, and the time weight is determined as W2 = exp(-Δt) 2 / (2σt 2 Let Δt be the time difference and σt = 30 days. This allows us to determine the atmospheric phase delay Φatm = Σ[WsWt(aΔP / Δh + bΔT)]. The coefficients a and b are solved using the least squares method, with an iteration threshold of 0.1 rad. Finally, the atmospheric phase is corrected by phase filtering using a modified Squeeze-and-Stretch filter to preserve the deformed phase. A quadratic polynomial fitting is then performed on the corrected phase to eliminate residual trend terms.

[0079] like Figure 4 As shown, in a preferred embodiment of the present invention, the step of performing dynamic deformation modeling specifically includes:

[0080] S301, extract the time series phase Δφ(t,r) from the atmospherically corrected interferogram set;

[0081] S302, perform linear deformation parameter inversion, and establish a system of linear equations for each pixel r: Δφl(ti,r)=4π / λ×v(r)×ti+∈i;

[0082] S303, subtract the linear deformation phase from the total phase to obtain the nonlinear deformation phase Δφn(t,r)=Δφ(t,r)-Δφl(t,r);

[0083] S304, perform piecewise nonlinear inversion, and establish a system of nonlinear equations for each pixel r: Δφn(ti,r)=4π / λ×A(r)×(1-e -ti / τ(r) )+∈i;

[0084] S305 generates a deformation rate cloud map and marks high-risk settlement areas; it also generates a deformation amplitude cloud map and a decay time constant cloud map.

[0085] In this embodiment of the invention, the total deformation phase is composed of the superposition of linear deformation phase and nonlinear deformation phase. The linear deformation phase reflects long-term slow deformation (such as subsidence in goaf areas) and changes linearly with time; the nonlinear deformation phase reflects instantaneous or short-term drastic deformation (such as blasting and landslides) and changes nonlinearly with time. First, the time series phase Δφ(t,r) is extracted from the atmospherically corrected interferometric atlas. Then, the linear deformation parameter is inverted. For each pixel r, a linear equation system is established: Δφl(ti,r)=4π / λ×v(r)×ti+∈i, where λ is the radar wavelength (e.g., 5.6cm), v(r) is the linear deformation rate, ∈i represents the residual term, and v(r) is solved by the least squares method. Then, subtracting the linear deformation phase from the total phase yields the nonlinear deformation phase Δφn(t,r) = Δφ(t,r) - Δφl(t,r), which allows for piecewise nonlinear inversion. For each pixel r, a system of nonlinear equations is established: Δφn(ti,r) = 4π / λ × A(r) × (1 - e -ti / τ(r) Let A(r) represent the amplitude A(r), estimated by nonlinear deformation phase peak value, and τ(r) represent the decay time constant. The initial range (e.g., 10–30 days) is set based on mining blasting records or landslide cases. Finally, a deformation rate cloud map is generated based on the linear deformation results to mark high-risk subsidence areas (e.g., deformation rate > 50 mm / yr); a deformation amplitude cloud map and a decay time constant cloud map are also generated based on the nonlinear deformation results to mark the temporal and spatial distribution of instantaneous deformation events (e.g., blasting, landslide).

[0086] like Figure 5 As shown, in a preferred embodiment of the present invention, the step of calculating the horizontal and vertical deformation components specifically includes:

[0087] S401 extracts satellite orbital parameters from the SAR image header file. The satellite orbital parameters include the ascending / descending orbit mode, incident angle, and azimuth angle.

[0088] S402, based on the digital elevation model, calculates the slope and aspect of each pixel. The slope represents the degree of inclination of the ground surface, and the aspect represents the direction of inclination of the ground surface.

[0089] S403 performs line-of-sight deformation projection decomposition. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle, while the vertical deformation component is determined by the incident angle.

[0090] In this embodiment of the invention, it should be noted that the deformation measured by InSAR is the projection value along the radar line-of-sight (LOS), which includes the combined contribution of vertical and horizontal deformation. The direction of the LOS deformation is determined by satellite orbital parameters (such as ascending or descending orbit) and the incident angle. The digital elevation model (DEM) provides elevation information of the mining area, which is used to calculate the terrain slope and aspect, and to assist in the decomposition of horizontal deformation. When determining the deformation components, satellite orbital parameters are first extracted from the SAR image header file. The satellite orbital parameters include the ascending / descending orbit mode, the incident angle, and the azimuth angle. Then, based on the digital elevation model, the slope and aspect of each pixel are calculated. The slope represents the degree of surface tilt, which affects the projection of horizontal deformation in the LOS direction. The aspect represents the direction of surface tilt, which determines the sign of the horizontal deformation projected onto the LOS. Then, the line-of-sight deformation projection decomposition is performed. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle. The sign of the projected deformation is determined by the relative relationship between the slope and the satellite flight direction (azimuth angle). The vertical deformation component is determined by the incident angle. The larger the vertical deformation, the larger the LOS deformation (positive correlation).

[0091] like Figure 6 As shown, this embodiment of the invention also provides a surface deformation monitoring system based on time-series InSAR, the system comprising:

[0092] The image screening module 100 is used to screen SAR images and generate interferograms, collect corresponding meteorological data and thermal infrared data, and retrieve the surface temperature through thermal infrared data.

[0093] The atmospheric phase correction module 200 is used to perform spatiotemporal adaptive atmospheric phase correction, taking the vertical pressure gradient and temperature anomaly as the driving factors of atmospheric delay, and separating the atmospheric phase through a spatiotemporal weighted model.

[0094] The deformation mechanism determination module 300 is used for dynamic deformation modeling, decomposing deformation into linear deformation and nonlinear deformation, and capturing the deformation mechanism of the mining area through segmented modeling.

[0095] The deformation component determination module 400 is used to retrieve the road network in the mining area as a geometric constraint, adopt the minimum cost flow algorithm, unwrap the coherent region, convert the unwrap phase into line-of-sight deformation, and combine InSAR line-of-sight deformation with digital elevation model to calculate the horizontal and vertical deformation components.

[0096] In a preferred embodiment of the present invention, the image screening module 100 includes:

[0097] The SAR image filtering unit is used to filter out SAR images whose spatial and temporal baselines match, and to remove SAR images with cloud cover greater than 10%.

[0098] Interferometric network building unit, used to construct temporal interferometric networks using a small baseline set strategy, and to perform multi-view processing on interferograms;

[0099] The coherence coefficient calculation unit is used to calculate the coherence coefficient threshold and remove low-coherence regions.

[0100] In a preferred embodiment of the present invention, the atmospheric phase correction module 200 includes:

[0101] The vertical pressure gradient unit is used to retrieve ERA5 data, extract the tropopause pressure PT and the surface pressure PS from the ERA5 data, and calculate the vertical pressure gradient ΔP / Δh=(PS-PT) / H, where H is the atmospheric elevation.

[0102] The temperature anomaly unit is used to calculate the residual between the surface temperature and the ERA5 temperature, identify temperature anomaly areas, and determine the temperature anomaly value ΔT.

[0103] Spatiotemporal weighted model unit, used to construct spatiotemporal weighted model and determine spatial weight W1 = exp(-d 2 / (2σs 2 )), d is the pixel spacing, σs = 500m, and the time weight is determined as W2 = exp(-Δt) 2 / (2σt 2 )), Δt is the time difference, σt = 30 days; determine the atmospheric phase delay Φatm = Σ[WsWt(aΔP / Δh+bΔT)], and solve for the coefficients a and b using the least squares method;

[0104] The atmospheric phase correction unit is used to correct the atmospheric phase, perform phase filtering, and perform quadratic polynomial fitting on the corrected phase to eliminate residual trend terms.

[0105] In a preferred embodiment of the present invention, the deformation mechanism determination module 300 includes:

[0106] The time series phase extraction unit is used to extract the time series phase Δφ(t,r) from the atmospherically corrected interferogram set.

[0107] The linear deformation inversion unit is used to perform linear deformation parameter inversion. For each pixel r, a system of linear equations is established: Δφl(ti,r)=4π / λ×v(r)×ti+∈i, where λ is the radar wavelength, v(r) is the linear deformation rate, and ∈i represents the residual term;

[0108] Nonlinear deformation phase unit, used to subtract linear deformation phase from total phase to obtain nonlinear deformation phase Δφn(t,r)=Δφ(t,r)-Δφl(t,r);

[0109] The nonlinear inversion unit is used to perform piecewise nonlinear inversion. For each pixel r, a system of nonlinear equations is established: Δφn(ti,r)=4π / λ×A(r)×(1-e -ti / τ(r) )+∈i, A(r) represents the amplitude A(r), and τ(r) represents the decay time constant;

[0110] The cloud map generation unit is used to generate deformation rate cloud maps and mark high-risk settlement areas; it also generates deformation amplitude cloud maps and decay time constant cloud maps.

[0111] In a preferred embodiment of the present invention, the deformation component determination module 400 includes:

[0112] The orbit parameter extraction unit is used to extract satellite orbit parameters from the SAR image header file. The satellite orbit parameters include the ascending / descending orbit mode, incident angle, and azimuth angle.

[0113] The slope and aspect calculation unit is used to calculate the slope and aspect of each pixel based on the digital elevation model. The slope represents the degree of inclination of the ground surface, and the aspect represents the direction of inclination of the ground surface.

[0114] The horizontal and vertical deformation elements are used to perform line-of-sight deformation projection decomposition. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle, while the vertical deformation component is determined by the incident angle.

[0115] The above description only details the preferred embodiments of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0116] It should be understood that although the steps in the flowcharts of the various embodiments of the present invention are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the various embodiments may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps.

[0117] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments described above. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and RAMbus dynamic RAM (RDRAM), etc.

[0118] Other embodiments of this disclosure will readily occur to those skilled in the art upon consideration of the disclosure in the specification and embodiments. This application is intended to cover any variations, uses, or adaptations of this disclosure that follow the general principles of this disclosure and include common knowledge or customary techniques in the art not disclosed herein. The specification and embodiments are to be considered exemplary only, and the true scope and spirit of this disclosure are indicated by the claims.

Claims

1. A method for monitoring surface deformation based on time-series InSAR, characterized in that, The method includes the following steps: SAR image screening and interferogram generation are performed, and corresponding meteorological and thermal infrared data are collected. The surface temperature is retrieved through thermal infrared data. Spatiotemporal adaptive atmospheric phase correction is performed, taking the vertical pressure gradient and temperature anomaly as the driving factors of atmospheric delay, and separating the atmospheric phase through a spatiotemporal weighted model; Dynamic deformation modeling is performed, decomposing deformation into linear and nonlinear deformation, and capturing the deformation mechanism of the mining area through segmented modeling; The mining area road network is used as a geometric constraint. The minimum cost flow algorithm is used to unwrap the coherent region and convert the unwrapped phase into line-of-sight deformation. The horizontal and vertical deformation components are calculated by combining InSAR line-of-sight deformation with digital elevation model.

2. The surface deformation monitoring method based on time-series InSAR according to claim 1, characterized in that, The steps for SAR image screening and interferogram generation specifically include: SAR images with matching spatial and temporal baselines were selected, and SAR images with cloud cover exceeding 10% were removed. A small baseline set strategy is adopted to construct a temporal interferometric network and perform multi-view processing on the interferogram; Calculate the coherence coefficient threshold and remove low-coherence regions.

3. The surface deformation monitoring method based on time-series InSAR according to claim 1, characterized in that, The steps for performing spatiotemporal adaptive atmospheric phase correction specifically include: Retrieve ERA5 data, extract PT (tropopause pressure) and PS (surface pressure) from ERA5 data, and calculate the vertical pressure gradient ΔP / Δh = (PS-PT) / H, where H is the atmospheric elevation. Calculate the residual between surface temperature and ERA5 temperature, identify temperature anomaly areas, and determine the temperature anomaly value ΔT; Construct a spatiotemporal weighted model and determine the spatial weight W1 = exp(-d 2 / (2σs 2 )), d is the pixel spacing, σs = 500m, and the time weight is determined as W2 = exp(-Δt) 2 / (2σt 2 )), Δt is the time difference, σt = 30 days; determine the atmospheric phase delay Φatm = Σ[WsWt(aΔP / Δh+bΔT)], and solve for the coefficients a and b using the least squares method; The atmospheric phase is corrected, phase filtering is performed, and the corrected phase is fitted with a quadratic polynomial to eliminate residual trend terms.

4. The surface deformation monitoring method based on time-series InSAR according to claim 1, characterized in that, The steps for performing dynamic deformation modeling specifically include: Extract the time series phase Δφ(t,r) from the atmospherically corrected interferogram set; To perform linear deformation parameter inversion, for each pixel r, establish a system of linear equations: Δφl(ti,r)=4π / λ×v(r)×ti+∈i, where λ is the radar wavelength, v(r) is the linear deformation rate, and ∈i represents the residual term; Subtracting the linear deformation phase from the total phase yields the nonlinear deformation phase Δφn(t,r) = Δφ(t,r) - Δφl(t,r); Perform piecewise nonlinear inversion, and establish a system of nonlinear equations for each pixel r: Δφn(ti,r)=4π / λ×A(r)×(1-e -ti / τ(r) )+∈i, A(r) represents the amplitude A(r), and τ(r) represents the decay time constant; Generate deformation rate cloud maps and mark high-risk settlement areas; generate deformation amplitude cloud maps and decay time constant cloud maps.

5. The surface deformation monitoring method based on time-series InSAR according to claim 1, characterized in that, The steps for calculating the horizontal and vertical deformation components specifically include: Extract satellite orbital parameters from the SAR image header file. These parameters include the ascending / descending orbit mode, incident angle, and azimuth angle. Based on the digital elevation model, the slope and aspect of each pixel are calculated. The slope represents the degree of inclination of the ground surface, and the aspect represents the direction of inclination of the ground surface. Line-of-sight deformation projection decomposition is performed. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle, while the vertical deformation component is determined by the incident angle.

6. A surface deformation monitoring system based on time-series InSAR, characterized in that, The system includes: The image filtering module is used to filter SAR images and generate interferograms, collect corresponding meteorological data and thermal infrared data, and retrieve the surface temperature through thermal infrared data. The atmospheric phase correction module is used for spatiotemporal adaptive atmospheric phase correction. It takes the vertical pressure gradient and temperature anomaly as the driving factors of atmospheric delay and separates the atmospheric phase through a spatiotemporal weighted model. The deformation mechanism determination module is used for dynamic deformation modeling, decomposing deformation into linear and nonlinear deformation, and capturing the deformation mechanism of the mining area through segmented modeling. The deformation component determination module is used to retrieve the road network in the mining area as a geometric constraint. It adopts the minimum cost flow algorithm to unwrap the coherent region, convert the unwrapped phase into line-of-sight deformation, and calculate the horizontal and vertical deformation components by combining the InSAR line-of-sight deformation with the digital elevation model.

7. The surface deformation monitoring system based on time-series InSAR according to claim 6, characterized in that, The image filtering module includes: The SAR image filtering unit is used to filter out SAR images whose spatial and temporal baselines match, and to remove SAR images with cloud cover greater than 10%. Interferometric network building unit, used to construct temporal interferometric networks using a small baseline set strategy, and to perform multi-view processing on interferograms; The coherence coefficient calculation unit is used to calculate the coherence coefficient threshold and remove low-coherence regions.

8. The surface deformation monitoring system based on time-series InSAR according to claim 6, characterized in that, The atmospheric phase correction module includes: The vertical pressure gradient unit is used to retrieve ERA5 data, extract the tropopause pressure PT and the surface pressure PS from the ERA5 data, and calculate the vertical pressure gradient ΔP / Δh=(PS-PT) / H, where H is the atmospheric elevation. The temperature anomaly unit is used to calculate the residual between the surface temperature and the ERA5 temperature, identify temperature anomaly areas, and determine the temperature anomaly value ΔT. Spatiotemporal weighted model unit, used to construct spatiotemporal weighted model and determine spatial weight W1 = exp(-d 2 / (2σs 2 )), d is the pixel spacing, σs = 500m, and the time weight is determined as W2 = exp(-Δt) 2 / (2σt 2 )), Δt is the time difference, σt = 30 days; determine the atmospheric phase delay Φatm = Σ[WsWt(aΔP / Δh+bΔT)], and solve for the coefficients a and b using the least squares method; The atmospheric phase correction unit is used to correct the atmospheric phase, perform phase filtering, and perform quadratic polynomial fitting on the corrected phase to eliminate residual trend terms.

9. The surface deformation monitoring system based on time-series InSAR according to claim 6, characterized in that, The deformation mechanism determination module includes: The time series phase extraction unit is used to extract the time series phase Δφ(t,r) from the atmospherically corrected interferogram set. The linear deformation inversion unit is used to perform linear deformation parameter inversion. For each pixel r, a system of linear equations is established: Δφl(ti,r)=4π / λ×v(r)×ti+∈i, where λ is the radar wavelength, v(r) is the linear deformation rate, and ∈i represents the residual term; Nonlinear deformation phase unit, used to subtract linear deformation phase from total phase to obtain nonlinear deformation phase Δφn(t,r)=Δφ(t,r)-Δφl(t,r); The nonlinear inversion unit is used to perform piecewise nonlinear inversion. For each pixel r, a system of nonlinear equations is established: Δφn(ti,r)=4π / λ×A(r)×(1-e -ti / τ(r) )+∈i, A(r) represents the amplitude A(r), and τ(r) represents the decay time constant; The cloud map generation unit is used to generate deformation rate cloud maps and mark high-risk settlement areas; it also generates deformation amplitude cloud maps and decay time constant cloud maps.

10. The surface deformation monitoring system based on time-series InSAR according to claim 6, characterized in that, The deformation component determination module includes: The orbit parameter extraction unit is used to extract satellite orbit parameters from the SAR image header file. The satellite orbit parameters include the ascending / descending orbit mode, incident angle, and azimuth angle. The slope and aspect calculation unit is used to calculate the slope and aspect of each pixel based on the digital elevation model. The slope represents the degree of inclination of the ground surface, and the aspect represents the direction of inclination of the ground surface. The horizontal and vertical deformation elements are used to perform line-of-sight deformation projection decomposition. The projection coefficient of the horizontal deformation component in the radar line-of-sight direction is determined by the slope and azimuth angle, while the vertical deformation component is determined by the incident angle.

Images