An engineering deformation and geological disaster monitoring method and system
By using a complementary fusion method of X-band and C-band InSAR data, the problems of ground hardware dependence and incomplete coverage of single bands in traditional monitoring technologies have been solved. This has enabled high-precision deformation monitoring of the entire area around large-scale projects, reduced monitoring costs in complex terrains, and improved coverage and early warning accuracy.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- CHINA UNIV OF GEOSCIENCES (WUHAN)
- Filing Date
- 2026-05-19
- Publication Date
- 2026-07-31
AI Technical Summary
Existing monitoring technologies rely on ground-based hardware facilities, and single-band coverage is incomplete, making it difficult to achieve high-precision deformation monitoring of the entire area surrounding large-scale projects. In particular, signals are prone to loss of coherence under complex terrain and vegetation cover, and there is a lack of unified benchmarks and atmospheric correction methods.
A complementary fusion method of X-band and C-band InSAR data is adopted. By acquiring multi-source remote sensing data, dual-band independent interferometric calculation is performed. Combined with precise orbital ephemeris and external digital elevation model, geocoding and atmospheric correction are performed to construct a fused deformation field. Atmospheric delay phase is removed by using spatiotemporal filtering algorithm to achieve high-precision monitoring without ground hardware support.
It has achieved high-precision, full-area, blind-spot-free monitoring in heterogeneous environments, reduced monitoring costs, ensured the spatial benchmark consistency of data fusion and the objectivity of early warning results, and increased the coverage from 15% to 92%.
Smart Images
Figure CN122218688B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of radar remote sensing mapping and engineering geological disaster monitoring technology, specifically relating to a method and system for monitoring engineering deformation and geological disasters. Background Technology
[0002] Safety monitoring of major projects such as large hydropower stations, cross-sea bridges, and deep-buried tunnels faces severe challenges. Monitoring needs have expanded from the single "engineering structure itself" to "the surrounding geological environment (such as reservoir bank landslides and access road slopes)".
[0003] Current deformation monitoring technologies have the following limitations: Traditional point-based monitoring is costly and difficult to maintain: Traditional GNSS (Global Navigation Satellite System) or total station monitoring requires the deployment of a large number of ground devices. In reservoir areas with steep terrain and dense vegetation or high mountain valleys, equipment installation is difficult and it is extremely susceptible to damage from landslides, power outages, or signal blockages, resulting in extremely high maintenance costs.
[0004] Single-band InSAR monitoring has physical blind spots: X-band (such as TerraSAR-X): Short wavelength (about 3cm), extremely high resolution, and very sensitive to tiny cracks and displacements on the dam surface. However, its penetration ability is weak, and once it encounters shrubs, grasslands or woodlands around the project, the signal is very prone to decorrelation, making potential landslide hazards in the surrounding area "invisible".
[0005] C-band (such as Sentinel-1): With a longer wavelength (approximately 5.6 cm), it has strong penetrating power and can obtain ground information through vegetation canopy, making it suitable for monitoring landslides. However, its resolution is relatively low (usually 10-20 m), which cannot meet the millimeter-level monitoring requirements of delicate artificial structures such as dams and bridge piers.
[0006] Multi-source fusion lacks a unified benchmark and atmospheric correction methods: Existing multi-band fusion technologies typically rely on ground-based GNSS stations to provide absolute coordinate benchmarks and remove atmospheric delay errors. If ground-based GNSS is eliminated, how to achieve high-precision geocoding alignment using only satellite data, and how to utilize algorithms to remove atmospheric phase interference from complex mountainous areas, are key technical challenges in achieving high-precision deformation monitoring for major engineering projects.
[0007] Therefore, there is an urgent need for a monitoring method that can combine the advantages of X-band and C-band and achieve atmospheric correction and high-precision fusion without the need for ground hardware support. Summary of the Invention
[0008] The purpose of this invention is to solve the problems of existing monitoring technologies relying on ground hardware facilities and incomplete coverage of a single band, and to provide a method and system for monitoring deformation of major engineering projects based on the complementary fusion of X-band and C-band InSAR data.
[0009] To achieve the above objectives, the present invention provides the following solution: A method for monitoring engineering deformation and geological hazards, comprising: Acquire and preprocess multi-source remote sensing data of the target monitoring area within the same monitoring period; wherein, the multi-source remote sensing data includes X-band single-view complex SAR image data, C-band single-view complex SAR image data, precise orbital ephemeris data, and external digital elevation model; Dual-band independent interferometry is performed on X-band single-view complex SAR image data and C-band single-view complex SAR image data to obtain dual-band solution results; wherein, the dual-band solution results include X-band line-of-sight deformation phase and C-band line-of-sight deformation phase; Based on precise orbital ephemeris data and external digital elevation models, the dual-band solution results are geocoded using a range-Doppler model, projected onto a grid in the same WGS-84 coordinate system, and atmospheric correction is performed on the dual-band solution results using a spatiotemporal filtering algorithm to obtain the clean deformation phase of the dual-band solution. Based on the dual-band pure deformation phase, the average coherence of each pixel in the grid in the X-band and C-band time series is calculated to construct the fused deformation field. By utilizing the statistical characteristics of the fused deformation field, a relative reference benchmark region is selected, and the deformation data of the entire field is relatively corrected. The deformation rate, cumulative deformation, and deformation acceleration of key monitoring points are extracted, and graded early warning information is output in combination with the preset engineering safety threshold.
[0010] Preferably, the method for obtaining the X-band line-of-sight deformation phase includes: The image with the largest comprehensive correlation coefficient in the X-band single-view complex SAR image data is selected as the main image, and the amplitude deviation index of all pixels in the main image is calculated. Pixels whose amplitude deviation index is less than a preset threshold are used as candidate points for permanent scatterers to construct a Delaunay triangulation. Linear regression analysis is performed on the phase difference between adjacent points in the Delaunay triangulation to obtain the X-band line-of-sight deformation phase, thus completing the solution of the X-band single-view complex SAR image data; wherein, the X-band line-of-sight deformation phase includes elevation error and linear deformation rate.
[0011] Preferably, the method for obtaining the C-band line-of-sight deformation phase includes: Based on the C-band single-view complex SAR image data, interferogram pairs are generated by using preset time baseline thresholds and spatial vertical baseline thresholds; The phase unwrapping of the interferogram pair is performed using the minimum cost flow algorithm to obtain the unwrapped phase; Based on the unwrapped phase, the deformation time series is inverted using the singular value decomposition method to obtain the optimal solution for the deformation rate; The optimal solution for the deformation rate is integrated over time to obtain the cumulative deformation at each time point, and a C-band line-of-sight deformation phase is generated.
[0012] Preferably, the method for atmospheric correction of the dual-band solution results using a spatiotemporal filtering algorithm includes: Based on the external digital elevation model, a linear regression model of unwrapped phase and elevation is constructed to estimate the vertical atmospheric phase related to elevation. The estimated vertical atmospheric phase is then subtracted from the dual-band solution results to obtain the residual phase after removing terrain-related errors. Using a preset time window, the residual phase is subjected to high-pass time filtering to separate the high-frequency time-domain components caused by atmospheric turbulence; The high-frequency components in the time domain are subjected to low-pass spatial filtering with a preset time window to extract the turbulent atmospheric phase screen that exhibits long-wavelength characteristics in space. The vertically stratified atmospheric phase and the turbulent atmospheric phase screen are simultaneously subtracted from the dual-band solution results to obtain the dual-band pure deformation phase.
[0013] Preferably, the method for constructing the fused deformation field includes: Based on the average coherence of time series, an adaptive fusion model containing three levels of logic, namely "structure priority", "blind spot filling" and "weighted transition", is constructed to generate an initial fusion deformation field. Spatially smooth the fusion boundary of the initial fusion deformation field to generate a seamless global surface deformation map, thus completing the construction of the fusion deformation field.
[0014] Preferably, the determination of the three-level logic includes: Structure priority logic: When the X-band coherence of a pixel is greater than or equal to the preset high confidence threshold, it is determined to be a structure-dominant region, and the weight coefficient of the X-band inverted deformation value is set to 1, and the weight coefficient of the C-band inverted deformation value is set to 0. Blind spot filling logic: When the X-band coherence of a pixel is less than a preset effective threshold and the C-band coherence is greater than or equal to the effective threshold, it is determined to be a monitoring blind spot, and the weight coefficient of the X-band inverted deformation value is set to 0, and the weight coefficient of the C-band inverted deformation value is set to 1. Weighted transition logic: When neither the structure priority logic nor the blind zone filling logic is satisfied, it is determined to be a transition zone, and the weights are allocated according to coherence normalization.
[0015] Preferably, the indicators for graded early warning include: rate indicators, incremental indicators, and tangent indicators.
[0016] The present invention also provides an engineering deformation and geological disaster monitoring system for implementing the method, comprising: The data access module is used to acquire and preprocess multi-source remote sensing data of the target monitoring area within the same monitoring period; wherein, the multi-source remote sensing data includes X-band single-view complex SAR image data, C-band single-view complex SAR image data, precise orbital ephemeris data, and external digital elevation models; The dual-frequency interferometric solution engine is used to perform independent interferometric solutions on X-band single-view complex SAR image data and C-band single-view complex SAR image data to obtain dual-frequency solution results; wherein, the dual-frequency solution results include X-band line-of-sight deformation phase and C-band line-of-sight deformation phase; The atmospheric correction module is used to geocode the dual-band solution results based on precise orbital ephemeris data and external digital elevation models using a range-Doppler model, project them onto a grid in the same WGS-84 coordinate system, and use a spatiotemporal filtering algorithm to perform atmospheric correction on the dual-band solution results to obtain clean deformation phases in the dual-band. The fusion module is used to calculate the average coherence of each pixel in the grid in the X-band and C-band time series based on the dual-band pure deformation phase, and to construct the fused deformation field. The early warning terminal is used to select a relative reference area using the statistical characteristics of the fused deformation field, perform relative correction on the deformation data of the whole field, extract the deformation rate, cumulative deformation and deformation acceleration of key monitoring points, and output graded early warning information in combination with preset engineering safety thresholds.
[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Achieved high-precision full-domain blind-spot-free monitoring in heterogeneous environments: Unlike traditional single-band monitoring or simple weighted average fusion, this invention adopts a three-level adaptive fusion model that includes "structure priority, blind spot filling, and weighted transition".
[0018] By setting a coherence threshold (0.75 / 0.25) for forced logical judgment, millimeter-level texture details of the X-band are forcibly preserved in the artificial structure area, avoiding noise pollution of low-resolution data; In vegetated areas, the penetration characteristics of the C-band were used to fill the decoherence blind zone of the X-band. This truly achieves seamless "point-to-surface" integrated monitoring of the dam body (high-precision points) and the reservoir bank slope (wide coverage).
[0019] 2. Eliminates dependence on ground hardware and significantly reduces monitoring costs in high mountain and canyon areas: To address the problem of difficulty in deploying GNSS and corner reflectors in complex terrain, this invention utilizes a spatiotemporal filtering algorithm based on a combination of high-pass time and low-pass space to mathematically separate and remove the atmospheric delay phase (APS).
[0020] This technology directly replaces the expensive ground-based hardware calibration system used in traditional monitoring. This frees the monitoring system from physical limitations such as power outages, equipment damage, or inaccessibility, greatly reducing the construction and maintenance costs of major projects throughout their entire lifecycle.
[0021] 3. Ensures the spatial reference consistency and objectivity of heterogeneous data fusion: This invention abandons the traditional registration method of manually selecting control points and adopts a geocoding technology that combines precise orbital ephemeris (POD) with external DEM.
[0022] This ensures that X-band and C-band data are strictly aligned at the pixel level in the WGS-84 coordinate system, eliminating geometric errors introduced by manual operation; The automated weight allocation formula based on coherence statistical characteristics avoids subjective interference from human experience in the fusion results, ensuring the objectivity and repeatability of the early warning results. Attached Figure Description
[0023] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are 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.
[0024] Figure 1 This is a flowchart illustrating the overall technical process of the present invention. Figure 2 This is a schematic diagram illustrating the principle of atmospheric correction; where, Figure 2 (a) in the diagram is the original interferometric phase diagram (including atmospheric errors); Figure 2 (b) in the diagram is a schematic diagram of spatiotemporal filtering separation (APS); Figure 2 (c) in the figure is the corrected deformation phase (pure) diagram; Figure 3 A schematic diagram illustrating the spatial complementarity and fusion effect of X-band and C-band; Figure 3 (a) in the image represents the X-band; Figure 3 (b) in the diagram represents the C-band; Figure 3 (c) in the diagram is a schematic of the fusion result; Figure 4 This is a logical architecture diagram of the monitoring system of the present invention; Figure 5 This is a logic diagram for tiered early warning decision-making. Detailed Implementation
[0025] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0026] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.
[0027] It should be noted that this embodiment uses a large-scale hydropower project (including a high arch dam and reservoir slope) located in a high mountain canyon area as an example for illustration, but the present invention is also applicable to deformation monitoring of major projects such as the slope of the entrance of a deep-buried tunnel, a cross-sea bridge, and a tailings dam.
[0028] Example 1: This embodiment adopts a fully software-based process, without relying on GNSS receivers or corner reflectors deployed on-site, and achieves monitoring entirely through closed-loop satellite remote sensing data.
[0029] Hardware environment: No on-site hardware required. The data processing unit uses a high-performance workstation (configuration: Intel Xeon processor, 128GB memory, NVIDIA RTX A6000 graphics card) to handle the heavy workload of SAR image interferometry processing and spatiotemporal filtering calculations.
[0030] Software environment: Based on the Linux operating system, integrating GAMMA or ISCE (InSAR Scientific Computing Environment) as the underlying processing core, and using Python to write automated fusion and batch processing scripts.
[0031] like Figure 1 As shown, a method for monitoring engineering deformation and geological hazards includes: S1: Acquire and preprocess multi-source remote sensing data of the target monitoring area within the same monitoring period; wherein, the multi-source remote sensing data includes X-band single-view complex (SLC) SAR image data, C-band single-view complex SAR image data, precise orbital ephemeris data (POD), and external digital elevation model (DEM). The preprocessing described in this embodiment includes registration and multi-view processing of the original images.
[0032] In this embodiment, the monitoring period is selected from January 1, 2024 to December 31, 2024.
[0033] X-band data (fine observations): Data source: TerraSAR-X satellite imagery. Mode: Spotlight mode, resolution 1m. 1m. Quantity: Acquire 25 images (revisited every 11 days), covering the main dam and power plant.
[0034] C-band data (wide-field observation): Data source: Sentinel-1A satellite imagery. Mode: Interferometric wide swath (IW) mode, resolution 5m. 20m. Quantity: Acquire 30 images (revisited every 12 days), covering a 20km range upstream and downstream of the reservoir area.
[0035] Unsupported positioning reference data: Precise orbit ephemeris (POD): Download the precise orbit files released by ESA and DLR, with a positioning accuracy better than 5cm, to replace ground control points (GCPs) for orbit correction. External terrain data: Download the SRTMV3 version of the 30m resolution digital elevation model (DEM) to remove terrain phase.
[0036] S2: Perform independent interferometry on X-band single-look complex SAR image data and C-band single-look complex SAR image data to obtain dual-band solution results; the dual-band solution results include X-band line-of-sight deformation phase and C-band line-of-sight deformation phase. Specifically, this step uses dual-thread parallel processing to extract deformation information of different ground features separately.
[0037] Thread A: X-band PS-InSAR processing (for artificial structures) includes: constructing a single master image interferometric pair for X-band images, employing Permanent Scatter Interferometry (PS-InSAR) technology, extracting high-coherence point targets using the amplitude deviation index, and calculating the line-of-sight deformation phase; a further implementation method for obtaining the X-band line-of-sight deformation phase includes: The image with the highest composite correlation coefficient among the X-band single-look complex SAR image data is selected as the master image. The specific calculation process is as follows: Assume the acquired X-band SAR image set contains... Scene images, each scene image ( Assuming the main image is assumed in turn, calculate its relationship with other images. The overall correlation coefficient of Jingcong images The formula for calculating the comprehensive correlation coefficient is: , In the formula, , , Images Compared to images Spatial vertical baseline, temporal baseline, and Doppler center frequency difference; , , These are the corresponding critical spatial baseline, critical temporal baseline, and critical Doppler center frequency difference, respectively. Calculate the values for all images. The image with the largest value is selected as the common master image (in this embodiment, an image from June 2024 is selected as the master image after calculation). Coarse registration is performed using precise orbital parameters, and fine registration is performed using intensity cross-correlation. The registration accuracy is controlled within 0.1 pixels.
[0038] Calculate the amplitude dispersion index for all pixels in the main image. The specific calculation method is as follows: for the registered... For a single-look complex SAR image, the amplitude stability of each pixel over time is statistically analyzed using the following formula: , , In the formula, For the first The amplitude value of that pixel in the image; This represents the average amplitude of the pixel over time. This represents the standard deviation of the time series amplitude of that pixel.
[0039] Pixels with amplitude deviation exponents less than a preset threshold are selected as candidate permanent scatterers to construct a Delaunay triangulation; specifically, a preset threshold is set. Highly stable permanent scattering bodies (PS) were screened on the dam surface and spillway concrete surface.
[0040] Linear regression analysis was performed on the phase difference between adjacent points in the Delaunay triangulation network to obtain the X-band line-of-sight deformation phase, thus completing the solution of X-band single-look complex SAR image data. The X-band line-of-sight deformation phase includes elevation error and linear deformation rate. Specifically, a phase observation equation was constructed. , Represents the total interference phase of the interferogram. This represents the phase contribution caused by surface deformation. This represents the phase contribution caused by topographic relief (i.e., topographic phase). This represents the phase contribution caused by atmospheric delay. The noise phase is represented by thermal noise and decoherent noise, etc. It is simulated using an SRTM DEM and subtracted. (Terrain phase).
[0041] Thread B: C-band SBAS-InSAR processing (for natural slopes) includes: constructing a multi-master image interferometry network for C-band imagery, employing Small Baseline Set Interferometry (SBAS-InSAR) technology to extract distributed targets, and calculating the line-of-sight deformation phase. A further implementation method for obtaining the C-band line-of-sight deformation phase includes: Based on C-band single-look complex SAR image data, interferogram pairs are generated by setting preset temporal baseline thresholds and spatial vertical baseline thresholds. In this embodiment, the temporal baseline threshold is set to 60 days, and the vertical spatial baseline threshold is set to 150 meters, generating approximately 80 pairs of interferograms, forming an interferometric network with high redundancy. A coherence threshold is also set. Even in the shrubland area along the reservoir bank, sufficient distributed scatterers (DS) can still be preserved by utilizing the penetrability of the C-band.
[0042] The minimum cost flow (MCF) algorithm is used to recover integer ambiguities, and the interferogram pairs are unwrapped to obtain the unwrapped phase. Specific steps include: Residual point identification: Calculate the phase gradient of the interferogram in the range and azimuth directions, and identify the phase discontinuities ("residues"), i.e., the points where the phase gradient is minimized by integrating along the smallest closed path. point.
[0043] Network Construction and Weight Setting: Construct a network connecting positive and negative residuals. Using the coherence graph generated in the previous steps as a weight reference, assign lower costs to paths with high coherence and higher costs to paths with low coherence. This step is to force phase transition "branch cuts" to preferentially pass through low-coherence areas (such as vegetation or water), thereby protecting the data quality of high-coherence areas.
[0044] Global optimization solution: Establish a global optimization model with the goal of minimizing total cost, solve for the optimal flow distribution in the network, and thus determine the integer blur degree of each pixel. .
[0045] Integral reconstruction: The determined integer number of cycles is added back into the wrapped phase, and the spatially continuous absolute unwrapped phase is obtained through path integration.
[0046] Based on the unwrapped phase, the deformation time series is inverted using the singular value decomposition (SVD) method to obtain the optimal solution for the deformation rate. The optimal solution for the deformation rate is then integrated over time to obtain the cumulative deformation at each time point, generating the C-band line-of-sight deformation phase. Specifically, the steps for inverting the deformation time series using SVD include: constructing an observation equation describing the functional relationship between phase and deformation rate. , in, for The design matrix is composed of time intervals between interferogram pairs; For demand 3D deformation rate vector; for dimensional unwrapped phase observation vector; design matrix for small baseline sets. The rank deficiency problem is solved by using singular value decomposition to calculate the matrix. generalized inverse matrix Thus, the optimal solution for the deformation rate is obtained: .
[0047] Finally, the deformation rate vector Perform time integration to obtain the cumulative deformation at each time point.
[0048] In the dual-band independent interferometric solution process, common steps include simulating the terrain phase using an external SRTM DEM and subtracting it from the interferometric phase to eliminate the influence of terrain on deformation inversion.
[0049] S3: Based on precise orbital ephemeris data and an external digital elevation model, the dual-band solution results are geocoded using a range-Doppler model and projected onto a raster grid in the same WGS-84 coordinate system. In this embodiment, a 10m grid covering the entire monitoring area is established. 10m standard grid.
[0050] The specific process includes: Construct a set of geometric positioning equations: For each pixel in the SAR image (corresponding to azimuth and time coordinates) and slant distance ), combined with satellite position vectors provided by precise orbits and velocity vector Establish the following system of simultaneous equations: Distance equation: Constraining the satellite and ground target points The geometric distance between them.
[0051] Doppler equations: ,in The center frequency of the Doppler wave. For radar wavelength, This represents the velocity of the ground target (usually set to 0).
[0052] Earth ellipsoid equation: Incorporating elevation data from an external DEM The target point is located on an ellipsoid defined by the Earth's semi-major and semi-minor axes.
[0053] Iterative solution for 3D coordinates: The Newton-Raphson iterative method is used to solve the above system of equations simultaneously to obtain the 3D rectangular coordinates of the target point in the Earth-centered Earth-fixed (ECEF) coordinate system. .
[0054] Coordinate transformation and resampling: The solved coordinates... Convert latitude and longitude to WGS-84 geodetic coordinate system The dual-band solution results are then projected onto a preset [projection method]. In a standard raster grid.
[0055] After geocoding and resampling, a spatiotemporal filtering algorithm is used to perform atmospheric correction on the dual-band solution results to obtain clean deformation phases for both bands. Specifically, the spatiotemporal filtering algorithm separates and removes atmospheric delay phases and residual orbital errors from the unwrapped phases of the two bands without relying on ground control points. Figure 2 As shown. Figure 2 (a) in the diagram is the original interferometric phase diagram (including atmospheric errors); Figure 2 (b) in the diagram is a schematic diagram of spatiotemporal filtering separation (APS); Figure 2 (c) in the figure is the corrected deformation phase (pure) diagram.
[0056] This embodiment addresses the significant atmospheric vertical stratification effect caused by large elevation differences in high mountain and canyon areas. Instead of the traditional single-filter method, it employs a two-stage correction strategy: "vertical stratification priority decoupling + turbulent hybrid residual filtering." A further implementation method for atmospheric correction of dual-band solution results using a spatiotemporal filtering algorithm includes: Level 1: Vertical Stratification Decoupling: Based on an external digital elevation model, a linear regression model of unwrapped phase and elevation is constructed to estimate the vertically stratified atmospheric phase related to elevation. The estimated vertically stratified atmospheric phase is then subtracted from the dual-band solution results to obtain the residual phase after removing terrain-related errors. Specifically, since the monitoring area is located in high mountains and valleys, atmospheric phase and elevation are strongly correlated. Before performing spatiotemporal filtering, the system first constructs a phase-elevation regression model: , in, To utilize pixel elevations obtained from an external DEM, This represents the residual phase. To prevent deformation points from interfering with atmospheric model estimations, this embodiment employs robust estimation or iterative weighted least squares to solve for the coefficients. and b.
[0057] in, denoted as , where is the phase-elevation correlation coefficient, physically representing the atmospheric phase delay gradient caused by a unit change in elevation; b is the intercept of the constant term, physically representing the atmospheric phase shift at the reference surface. To prevent deformation points (as statistical outliers) from interfering with atmospheric model estimations, this embodiment employs Iterative Reweighted Least Squares (IRLS) to solve for the coefficients. And b, the specific steps include: Step 1: Initial estimation. A linear fit is performed on all pixels using ordinary least squares (OLS) to obtain... And the initial estimate of b.
[0058] Step 2: Residual Calculation and Weight Assignment. Calculate the residual between the observed phase and the model-predicted phase for each pixel. The weights of each point are determined using IGGIII or Huber weight functions. The specific strategy is: when the residual When the value is less than a preset threshold (mainly caused by atmospheric factors), a high weight is assigned; when the residual... When the value exceeds the threshold (mainly caused by surface deformation), it is identified as an anomaly and assigned minimal or even zero weight.
[0059] Step 3: Weighted Iteration. Based on the weight matrix. Resolve the system of linear equations ,renew and The value; Step 4: Convergence Determination. Repeat steps 2 and 3 until the parameter difference between two consecutive iterations is less than the preset tolerance (e.g., ...). ), output the final and the value of b. Using the finally obtained coefficients... b) Estimate and separate and remove the topographically related vertical atmospheric phase from the original interferometric phase. Prioritize separating and removing the topographically related vertical atmospheric phase from the original interferometric phase.
[0060] Second stage: Turbulent Mixing Filtering: Using a preset time window, high-pass time filtering is applied to the residual phase to separate the high-frequency components in the time domain caused by atmospheric turbulence; low-pass spatial filtering is applied to the high-frequency components in the time domain within a preset time window to extract the turbulent atmospheric phase screen that exhibits long-wavelength characteristics in space; the vertically stratified atmospheric phase and the turbulent atmospheric phase screen are simultaneously subtracted from the dual-band solution results to obtain the dual-band pure deformable phase.
[0061] Specifically, regarding the residual phase after removing vertical stratification (Mainly caused by atmospheric turbulence), separation is achieved by utilizing its stochastic characteristics in the spatiotemporal domain: 1. High-pass time filtering: Set the time window to 365 days, perform high-pass filtering on the residual phase sequence to filter out non-atmospheric signals that change slowly with the seasons, and retain the high-frequency fluctuating atmospheric turbulence components.
[0062] 2. Low-pass spatial filtering: The above results are subjected to spatial domain low-pass filtering (window set to 1-2km) to extract the pure turbulent atmospheric phase screen (APS) with long wavelength correlation in space.
[0063] 3. Final correction: The vertical stratification phase estimated in the first stage and the turbulence phase extracted in the second stage are combined into the total atmospheric error and subtracted from the original unwrapped phase to obtain a pure deformation phase with a high signal-to-noise ratio.
[0064] S4: Based on the dual-band pure deformation phase, calculate the average coherence of each pixel in the raster grid in the X-band and C-band time series, and construct the fused deformation field. For example... Figure 3 As shown. Figure 3 (a) in the image represents the X-band; Figure 3 (b) in the diagram represents the C-band; Figure 3 (c) in the diagram is a schematic diagram of the fusion result.
[0065] A further implementation method involves constructing a fused deformation field, comprising: Based on the average coherence of time series data, an adaptive fusion model is constructed, incorporating three levels of logic: "structure priority," "blind spot filling," and "weighted transition," to generate an initial fused deformation field. The fusion boundary of the initial fused deformation field is spatially smoothed to generate a seamless global surface deformation map, thus completing the construction of the fused deformation field. For any pixel in the grid, the fused deformation rate is calculated using a linear weighting formula. : ; in, , These are the inversion deformation values for the X-band and C-band, respectively. , These are the weighting coefficients; In this implementation, the automatic weight distribution algorithm based on coherence characteristics is used to assign weights to each pixel within the grid. Perform a traversal and read the corresponding X-band coherence. Coherence with C-band The weighting coefficients are automatically calculated based on the preset judgment logic. , The specific algorithm flow is as follows: A further implementation method involves the following: the three-level logic determination includes: (1) Structure priority logic (corresponding to structure priority region): When the X-band coherence of a pixel is greater than or equal to the preset high confidence threshold, it is determined to be a structure-dominant region, and the weight coefficient of the X-band inverted deformation value is set to 1, and the weight coefficient of the C-band inverted deformation value is set to 0; Specifically, firstly, it is determined that... Is it greater than or equal to the high confidence threshold? (In this embodiment, it is set to 0.75). If the following conditions are met ( The point is determined to be either an artificial structure or a highly stable area of bare rock. To preserve high-resolution texture details in the X-band to the maximum extent, the algorithm forcibly locks the weights: .
[0066] (2) Blind spot filling logic (corresponding to effective supplementary area): When the X-band coherence of a pixel is less than the preset effective threshold and the C-band coherence is greater than or equal to the effective threshold, it is determined to be a monitoring blind spot. The weight coefficient of the X-band inverted deformation value is set to 0, and the weight coefficient of the C-band inverted deformation value is set to 1. Specifically, if step (1) is not satisfied, the system continues to determine whether the " and If the following conditions are met: the point is determined to be a vegetation cover area that is X-band decoherent but C-band valid. The algorithm then uses C-band data for blind spot filling. .
[0067] (3) Weighted Transition Logic (corresponding to Weighted Transition Zone): When neither the structure priority logic nor the blind zone filling logic is satisfied, it is determined to be a transition zone, and the weights are allocated according to coherence normalization. Specifically, if neither of the above two conditions is satisfied, the point is determined to be a mixed feature transition zone. The algorithm performs coherence-based normalized weighted calculation to ensure the smoothness of the fusion boundary: .
[0068] (4) Fusion rate generation: Substitute the calculated weights into the linear fusion formula. The final fusion deformation rate of the pixel is obtained.
[0069] S5: Utilize the statistical characteristics of the fused deformation field to select a relative reference area, perform relative correction on the full-field deformation data, and extract the deformation rate v, cumulative deformation S, and deformation acceleration at key monitoring points. It outputs graded early warning information based on preset engineering safety thresholds.
[0070] The selection criterion for the relative reference base area mentioned in step S5 is: using the fused full-time average coherence map, the following areas are selected. A high-coherence pixel set was selected; combined with external geological vector data, only pixels located in bedrock outcrops and more than 2 kilometers away from the engineering excavation face or reservoir bank were retained; the standard deviation of the deformation time series of the selected points was calculated. Remove Outliers identified by doubling the mean square error; the mean deformation of the remaining pixels is defined as the "zero deformation baseline," which is subtracted from the deformation sequence of all monitored points. For example... Figure 5 As shown.
[0071] A further implementation method involves calculating the following indicators for tiered early warning: Rate index ( ): The absolute value of the annual average deformation rate at the monitoring point; Acceleration index ( ): The acceleration of the deformation curve at the monitoring point is used to characterize the acceleration trend of deformation; Cumulative deformation index ( The cumulative deformation of the monitoring point within a specified period. Simultaneously, according to engineering design specifications, the safety threshold for the cumulative deformation of this monitoring point is set as follows: .
[0072] Based on the above indicators, the system's tiered early warning logic is as follows: 1) Red Alert: When judging cumulative deformation variables The system will prompt an immediate stop and issue an alarm when the deformation is in a state of continuous acceleration. 2) Orange warning: When the conditions for a red alert are not met, but the judgment rate is high. or acceleration When triggered, the system prompts that encrypted monitoring of the area is required; 3) Yellow alert: When the conditions for an orange warning are not met, but the judgment rate is... When triggered, the system prompts that data analysis needs to be strengthened and continuous monitoring should be maintained; 4) Green Normal: When none of the above conditions are met, the system determines that it is within the normal and safe range.
[0073] The selection of relative reference benchmark areas and early warning decisions include: 1. Selection of virtual reference points: The algorithm automatically searches for areas on the fused map that are more than 3 km away from the dam, marked as "granite bedrock outcrops" on the geological map, and have an InSAR coherence of more than 0.9 throughout the entire time period.
[0074] The average value of 100 pixels in the region is selected as the "zero deformation reference", and the deformation values of all monitoring points are corrected relative to this reference.
[0075] 2. Tiered early warning implementation: Monitoring results: The average annual relative deformation rate v of the middle part of the dam crest is 3.5 mm / yr (within the normal range of thermal expansion and contraction).
[0076] The relative deformation rate v was 42 mm / yr on the slope 1.5 km from the dam on the left bank of the reservoir.
[0077] Warning triggered: For dams: When the velocity v < 10 mm / yr, the system displays green (safe).
[0078] For the left bank slope: velocity v > 20 mm / yr, and acceleration over the past month An anomaly was detected, triggering an orange alert (Level II) on the system.
[0079] Decision support: The system automatically generates reports suggesting that the engineering management department conduct on-site inspections of the left bank slope, while the dam maintains routine inspections.
[0080] Through the above implementation method: Improved coverage: Compared to using the X-band alone, the monitoring coverage has increased from 15% (covering only buildings) to 92% (covering the entire watershed).
[0081] Accuracy verification: Although there is no GNSS for absolute correction, the standard deviation of the X-band on stable rock mass is controlled within ±2mm through atmospheric filtering and relative benchmark, which fully meets the needs of engineering survey and trend early warning.
[0082] Example 2: like Figure 4 As shown, the present invention also provides an engineering deformation and geological disaster monitoring system for implementing the method of embodiment one, comprising: The data access module is used to acquire and preprocess multi-source remote sensing data of the target monitoring area within the same monitoring period. The multi-source remote sensing data includes X-band single-view complex SAR image data, C-band single-view complex SAR image data, precise orbital ephemeris data, and external digital elevation models.
[0083] The dual-frequency interferometric solution engine is used to perform independent interferometric solutions on X-band single-view complex SAR image data and C-band single-view complex SAR image data to obtain dual-frequency solution results. The dual-frequency solution results include X-band line-of-sight deformation phase and C-band line-of-sight deformation phase.
[0084] The atmospheric correction module is used to geocode the dual-band solution results based on precise orbital ephemeris data and external digital elevation models using a distance-Doppler model, project them onto a raster grid in the same WGS-84 coordinate system, and use a spatiotemporal filtering algorithm to perform atmospheric correction on the dual-band solution results to obtain clean deformation phases for both frequencies.
[0085] The fusion module is used to calculate the average coherence of each pixel in the raster grid in the X-band and C-band time series based on the dual-band pure deformation phase, and to construct the fused deformation field. Figure 4 The fusion correction model in the model consists of the atmospheric correction module and the fusion module.
[0086] The early warning terminal is used to select a relative reference area by utilizing the statistical characteristics of the fused deformation field, perform relative correction on the deformation data of the entire field, extract the deformation rate, cumulative deformation and deformation acceleration of key monitoring points, and output graded early warning information in combination with preset engineering safety thresholds.
[0087] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.
Claims
1. An engineering deformation and geological disaster monitoring method, characterized in that, include: Acquire and preprocess multi-source remote sensing data of the target monitoring area within the same monitoring period; wherein, the multi-source remote sensing data includes X-band single-view complex SAR image data, C-band single-view complex SAR image data, precise orbital ephemeris data, and external digital elevation model; Dual-band independent interferometry is performed on X-band single-view complex SAR image data and C-band single-view complex SAR image data to obtain dual-band solution results; wherein, the dual-band solution results include X-band line-of-sight deformation phase and C-band line-of-sight deformation phase; Based on precise orbital ephemeris data and external digital elevation models, the dual-band solution results are geocoded using a range-Doppler model, projected onto a grid in the same WGS-84 coordinate system, and atmospheric correction is performed on the dual-band solution results using a spatiotemporal filtering algorithm to obtain the clean deformation phase of the dual-band solution. Based on the dual-band pure deformation phase, the average coherence of each pixel in the grid in the X-band and C-band time series is calculated to construct the fused deformation field. By utilizing the statistical characteristics of the fused deformation field, a relative reference benchmark region is selected, and the deformation data of the entire field is relatively corrected. The deformation rate, cumulative deformation and deformation acceleration of key monitoring points are extracted, and graded early warning information is output in combination with the preset engineering safety threshold. Methods for constructing fused deformation fields include: Based on the average coherence of time series, an adaptive fusion model containing three levels of logic, namely "structure priority", "blind spot filling" and "weighted transition", is constructed to generate the initial fusion deformation field; Spatially smooth the fusion boundary of the initial fused deformation field to generate a seamless global surface deformation map, thus completing the construction of the fused deformation field; The three-level logic determination includes: Structure priority logic: When the X-band coherence of a pixel is greater than or equal to the preset high confidence threshold, it is determined to be a structure-dominant region, and the weight coefficient of the X-band inverted deformation value is set to 1, and the weight coefficient of the C-band inverted deformation value is set to 0. Blind spot filling logic: When the X-band coherence of a pixel is less than a preset effective threshold and the C-band coherence is greater than or equal to the effective threshold, it is determined to be a monitoring blind spot, and the weight coefficient of the X-band inverted deformation value is set to 0, and the weight coefficient of the C-band inverted deformation value is set to 1. Weighted transition logic: When neither the structure priority logic nor the blind zone filling logic is satisfied, it is determined to be a transition zone, and the weights are allocated according to coherence normalization.
2. The method for monitoring engineering deformation and geological disasters according to claim 1, characterized in that, Methods for obtaining the X-band line-of-sight deformation phase include: The image with the largest comprehensive correlation coefficient in the X-band single-view complex SAR image data is selected as the main image, and the amplitude deviation index of all pixels in the main image is calculated. Pixels whose amplitude deviation index is less than a preset threshold are used as candidate points for permanent scatterers to construct a Delaunay triangulation. Linear regression analysis is performed on the phase difference between adjacent points in the Delaunay triangulation to obtain the X-band line-of-sight deformation phase, thus completing the solution of the X-band single-view complex SAR image data; wherein, the X-band line-of-sight deformation phase includes elevation error and linear deformation rate.
3. The method for monitoring engineering deformation and geological hazards according to claim 1, characterized in that, Methods for obtaining the C-band line-of-sight deformation phase include: Based on the C-band single-view complex SAR image data, interferogram pairs are generated by using preset time baseline thresholds and spatial vertical baseline thresholds; The phase unwrapping of the interferogram pair is performed using the minimum cost flow algorithm to obtain the unwrapped phase; Based on the unwrapped phase, the deformation time series is inverted using the singular value decomposition method to obtain the optimal solution for the deformation rate; The optimal solution for the deformation rate is integrated over time to obtain the cumulative deformation at each time point, and a C-band line-of-sight deformation phase is generated.
4. The method for monitoring engineering deformation and geological hazards according to claim 1, characterized in that, The method for atmospheric correction of the dual-band solution results using a spatiotemporal filtering algorithm includes: Based on the external digital elevation model, a linear regression model of unwrapped phase and elevation is constructed to estimate the vertical atmospheric phase related to elevation. The estimated vertical atmospheric phase is then subtracted from the dual-band solution results to obtain the residual phase after removing terrain-related errors. Using a preset time window, the residual phase is subjected to high-pass time filtering to separate the high-frequency time-domain components caused by atmospheric turbulence; The high-frequency components in the time domain are subjected to low-pass spatial filtering with a preset time window to extract the turbulent atmospheric phase screen that exhibits long-wavelength characteristics in space. The vertically stratified atmospheric phase and the turbulent atmospheric phase screen are simultaneously subtracted from the dual-band solution results to obtain the dual-band pure deformation phase.
5. The method for monitoring engineering deformation and geological hazards according to claim 1, characterized in that, The indicators for graded early warning include: rate indicators, acceleration indicators, and cumulative deformation indicators.
6. An engineering deformation and geological hazard monitoring system, used to implement the engineering deformation and geological hazard monitoring method according to any one of claims 1-5, characterized in that, include: The data access module is used to acquire and preprocess multi-source remote sensing data of the target monitoring area within the same monitoring period; wherein, the multi-source remote sensing data includes X-band single-view complex SAR image data, C-band single-view complex SAR image data, precise orbital ephemeris data, and external digital elevation models; The dual-frequency interferometric solution engine is used to perform independent interferometric solutions on X-band single-view complex SAR image data and C-band single-view complex SAR image data to obtain dual-frequency solution results; wherein, the dual-frequency solution results include X-band line-of-sight deformation phase and C-band line-of-sight deformation phase; The atmospheric correction module is used to geocode the dual-band solution results based on precise orbital ephemeris data and external digital elevation models using a range-Doppler model, project them onto a grid in the same WGS-84 coordinate system, and use a spatiotemporal filtering algorithm to perform atmospheric correction on the dual-band solution results to obtain clean deformation phases in the dual-band. The fusion module is used to calculate the average coherence of each pixel in the grid in the X-band and C-band time series based on the dual-band pure deformation phase, and to construct the fused deformation field. The early warning terminal is used to select a relative reference area using the statistical characteristics of the fused deformation field, perform relative correction on the deformation data of the whole field, extract the deformation rate, cumulative deformation and deformation acceleration of key monitoring points, and output graded early warning information in combination with preset engineering safety thresholds.