Surface deformation monitoring method based on multi-polarization time-series SAR data

Abstract

The present application relates to the technical field of remote sensing monitoring, in particular to a ground surface deformation monitoring method based on multi-polarization time-series SAR data, comprising: collecting multi-temporal multi-polarization SAR images of the same orbit in a monitoring area, and obtaining a time-series multi-polarization SAR data set through radiation calibration, filtering and registration preprocessing; extracting stable scatterer point targets by using a PS-InSAR algorithm, optimizing time-series phase unwrapping in combination with multi-polarization scattering characteristics, eliminating atmospheric delay, orbit and residual terrain phase error, and obtaining an optimized deformation phase sequence; separating linear deformation rate and nonlinear deformation time series through deformation modeling, and generating a continuous ground surface deformation field through spatial interpolation. The method expands the analysis dimension of radar data, effectively suppresses multiple types of phase interference, improves deformation inversion accuracy, and realizes fine and continuous monitoring of regional ground surface deformation full coverage.

Description

Technical Field

[0001] This invention relates to the field of remote sensing monitoring technology, and in particular to a method for monitoring surface deformation based on multipolar time-series SAR data. Background Technology

[0002] Surface deformation monitoring is widely used in geological environment exploration and land space management. Conventional large-scale surface deformation monitoring generally uses single-polarization time-series SAR data combined with PS-InSAR algorithms for calculation and analysis. Conventional data processing only uses single-polarization radar imagery for calculation. Core steps such as stable scatterer point target identification and time-series phase unwrapping are all completed based on single-dimensional scattering information. The overall data analysis dimension is limited, and deformation calculation is limited to the inversion of discrete point parameters.

[0003] Single-polarization SAR data can only capture limited surface scattering features. Under complex surface environments, interferometric phases are prone to superimposed atmospheric delay errors, orbital errors, and residual topographic phase errors. Traditional phase unwrapping methods lack multi-dimensional data correction capabilities and cannot simultaneously remove various interfering phase components. Time-series phase sequences are prone to retaining system biases, directly affecting the basic accuracy of subsequent deformation inversion.

[0004] Traditional deformation calculation architectures struggle to decompose and analyze differentiated deformation components, failing to independently distinguish the temporal variation patterns of linear and nonlinear deformations. Discretely distributed monitoring points cannot achieve full regional coverage, resulting in spatial gaps in the monitoring range. It is necessary to expand the polarization utilization of SAR data, improve phase unwrapping optimization mechanisms, decompose deformation parameters of different attributes, and address the spatial coverage shortcomings of discrete point monitoring to meet the needs for continuous and refined monitoring of surface deformation at regional scales. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by proposing a surface deformation monitoring method based on multipolar time-series SAR data.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a surface deformation monitoring method based on multi-polarization time-series SAR data, comprising: Collect multi-temporal and multi-polarization synthetic aperture radar (SAR) images of the preset monitoring area and the same orbit, and perform radiometric calibration, filtering and registration preprocessing on all SAR images to obtain a time-registered multi-polarization SAR dataset. Based on the time-registered multipolar SAR dataset, stable scatterer point targets within the monitoring area are identified and extracted using the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm. For the extracted stable scatterer point targets, the temporal phase unwrapping of the stable scatterer point targets is optimized by utilizing the multi-polarization scattering characteristics to obtain an optimized deformation phase sequence that removes atmospheric delay error, orbital error and residual terrain phase error; Geocoding and deformation modeling are performed on the optimized deformation phase sequence to separate and extract the surface linear deformation rate and nonlinear deformation time series of the monitoring area; Spatial interpolation is performed on the extracted linear and nonlinear deformation time series of the land surface to generate a continuously spatially distributed land surface deformation field in the monitoring area, and the deformation monitoring results are output.

[0007] As a further aspect of the present invention, the step of performing radiometric calibration, filtering, and registration preprocessing on all SAR images to obtain a time-registered multi-polarization SAR dataset specifically involves: All acquired multi-temporal, multi-polarization SAR images are sequentially radiometrically calibrated to convert the original echo signal intensity into backscattering coefficients, thereby eliminating sensor system errors. For the radiometrically calibrated SAR images of each polarization channel, an adaptive filtering algorithm is applied to suppress the speckle noise of the images while preserving the edge and texture information of the images. A common master image is selected as the SAR image with the best imaging quality in the middle of the time series, and all other secondary images in the time series are registered with the common master image. During registration, coarse registration is performed based on orbital parameters and the digital elevation model (DEM), and fine registration is performed at the pixel level using cross-correlation or feature matching methods to ensure that the registration error between all slave images and the common master image is less than a preset threshold. For image pairs that have completed pixel-level registration, the spectral diversity method is further used to optimize the sub-pixel-level registration and obtain the registration offset field. All images are resampled using the registration offset field to geometrically align SAR images of all time phases, forming a time-registered multi-polarization SAR dataset.

[0008] As a further aspect of the present invention, based on the time-registered multi-polarization SAR dataset, stable scatterer point targets within the monitoring area are identified and extracted using the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm. Specifically: For the time-registered multi-polarization SAR dataset, for each pixel location, the amplitude deviation index and phase stability index of its backscattering coefficient time series under all time phases and all polarization channels are calculated; Combining the amplitude deviation index and the phase stability index, a multidimensional feature space is constructed, and the position of the pixel in the multidimensional feature space is used as a measure of its scattering stability. Based on the multidimensional feature space, a clustering analysis algorithm is used to automatically identify pixels that simultaneously have low amplitude deviation index and high phase stability index as candidate permanent scattering points. For the identified candidate permanent scatterers, their dominant scattering mechanisms are distinguished by their scattering matrix information under different polarization channels. By combining the preset stability threshold of scattering mechanism, connectivity analysis is performed on candidate permanent scatterer points belonging to the same scattering mechanism, isolated noise points are eliminated, and finally a set of stable scatterer point targets that are spatially coherently distributed within the monitoring area is extracted.

[0009] As a further aspect of the present invention, the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm specifically includes: A differential interferometric network suitable for multipolar data is constructed. The differential interferometric network adds interferometric pairs that consider polarization basis changes on the basis of conventional short-spatial-temporal baseline combinations. For each differential interferogram in the differential interferometry network, on the stable scatterer point target, the phase deviation introduced by the polarization basis inconsistency is estimated and corrected using the multi-polarization information in the time-registered multi-polarization SAR dataset to form the polarization-corrected differential interferometric phase. The terrain phase contribution on the stable scatterer point target is simulated and removed using external digital elevation model (DEM) data to generate residual differential interferometric phase. For the residual differential interference phase, a graph-based three-dimensional phase unwrapping method is used to jointly unwrap the phase in the spatial and temporal domains to obtain the absolute differential interference phase; Based on the spatiotemporal correlation of atmospheric phase, atmospheric phase components are separated from the absolute differential interferometric phase using a filtering method. The residual digital elevation model error is estimated by using spatially continuous and temporally uncorrelated high-frequency phase residuals. The atmospheric phase component and the residual digital elevation model error are sequentially removed from the absolute differential interferometric phase to obtain the optimized residual differential interferometric phase. The optimized residual differential interferometric phase sequence is the optimized deformation phase sequence.

[0010] As a further aspect of the present invention, the step of optimizing the temporal phase unwrapping of the extracted stable scattering point target using multi-polarization scattering characteristics specifically involves: For each stable scattering point target, extract its scattering matrix or coherence matrix under all time phases and all polarization channels; The scattering matrix or coherence matrix is ​​decomposed into eigenvalues ​​to obtain the dominant scattering mechanism and corresponding scattering contribution of the stable scattering point target in different time phases. Based on the dominant scattering mechanism and the time-series stability of the corresponding scattering contribution, the reliability of the phase of the stable scattering point target in different interferometric pairs is evaluated, and a weighting coefficient based on polarization information is assigned to the phase of each interferometric pair. During the phase unwrapping process, the weighting coefficients based on polarization information are introduced into the smoothing constraint or cost function of unwrapping, so that the phase unwrapping process tends to follow the change path of high-weight, high-reliability phases. For the unwrapped absolute phase, additional constraints are constructed using the phase relationship between the multi-polarization channels. Cross-validation and consistency correction are performed on the unwrapping results to correct unwrapping errors, thereby optimizing the phase unwrapping results and reducing the uncertainty of phase unwrapping.

[0011] As a further aspect of the present invention, the optimized deformation phase sequence is geocoded and modeled for deformation, and the linear deformation rate and nonlinear deformation time series of the monitored area are separated and extracted, specifically as follows: Using the orbital parameters and imaging geometry of the radar system, the optimized deformation phase sequence is converted from the slant range geometry of the radar coordinate system to the ground geometry in the geographic coordinate system, thus completing the geocoding. Assuming that deformation consists of a linear trend term, a seasonal periodic term, and an irregular nonlinear residual term, a time series model of surface deformation is constructed. Using the geocoded optimized deformation phase sequence as observations, the parameters of the time series model of the surface deformation are solved by least squares estimation or singular value decomposition. The coefficients of the linear trend term are directly extracted from the solved model parameters and used as the linear deformation rate of the land surface. Subtracting the fitted linear trend term and the known seasonal periodic term from the optimized deformation phase sequence, the remaining part is the nonlinear deformation time series, which includes unmodeled deformation signals, residual noise, and abrupt deformation information.

[0012] As a further aspect of the present invention, the step of spatially interpolating the extracted linear and nonlinear deformation time series of the land surface to generate a continuously spatially distributed land surface deformation field in the monitoring area specifically involves: The geographical location of the extracted stable scatterer point target is used as the known sample point, and the corresponding surface linear deformation rate is used as the sample value. Based on the natural neighbor method, Kriging interpolation method or inverse distance weighting method, spatial interpolation calculation is performed on the surface linear deformation rate of the known sample points to generate a continuous spatial distribution of linear deformation rate raster map covering the entire monitoring area. For each time point in the nonlinear deformation time series, the nonlinear deformation value of the stable scattering point target at the corresponding time point is taken as a sample value; Using the same spatial interpolation method as the linear deformation rate, independent spatial interpolation is performed on the nonlinear deformation sample at each time point to generate a nonlinear deformation spatial distribution raster map corresponding to each time point. The continuous spatial distribution of linear deformation rate raster map is integrated with the nonlinear deformation spatial distribution raster map at all time points to form a complete surface deformation field. The surface deformation field is continuous in space and contains linear trends and nonlinear details in time.

[0013] As a further aspect of the present invention, the connectivity analysis of candidate permanent scatterer points belonging to the same scattering mechanism, combined with a preset scattering mechanism stability threshold, to eliminate isolated noise points, specifically involves: A stability threshold for the scattering mechanism is set over a time series, which defines the maximum allowable range of variation of the scattering mechanism between adjacent time phases; Calculate the time series of scattering mechanisms for each candidate permanent scatterer point and analyze the changes in scattering mechanisms between adjacent time phases; Candidate points whose changes are consistently below the stability threshold of the scattering mechanism are identified as permanent scattering points with a stable scattering mechanism. In the spatial domain, connectivity analysis of eight-neighbor or four-neighbor domains is performed on all permanent scattering points determined to have stable scattering mechanisms. Identify and mark all spatially connected point groups, and discard isolated point groups or single points with fewer than a preset minimum number of pixels as noise points. The group of connected points with a number of pixels greater than or equal to a preset lower limit is retained as the final extracted stable scattering point target.

[0014] As a further aspect of the present invention, during the phase unwrapping process, the weighting coefficients based on polarization information are introduced into the smoothing constraint or cost function of the unwrapping process, specifically as follows: Define a cost function for phase unwrapping, which includes the weighted sum of squares of the differences between the unwrapped phase gradient and the wrapped phase gradient; The weighting coefficients based on polarization information are used as weighting factors for the corresponding interference pairs and point targets in the cost function. The weighting factor for high reliability phases is large, and the weighting factor for low reliability phases is small. When performing phase unwrapping using the cost function minimization, high-weighted interference contributes more to the phase gradient difference in the total cost, forcing the unwrapping result to be more inclined to satisfy the continuity of the high-reliability phase. By using an iterative optimization algorithm, the unwrapped phase field that minimizes the total weighted cost is found, and the optimized phase unwrapping result is obtained.

[0015] As a further aspect of the present invention, the spatial interpolation calculation of the surface linear deformation rate of the known sample points based on the natural neighbor method, Kriging interpolation method, or inverse distance weighting method specifically includes: If the natural neighbor method is used, for each pixel to be interpolated in the monitoring area, its natural neighbor is found among the known sample points. The natural neighbor is the sample point in the Thiessen polygon that is adjacent to the polygon in which the pixel to be interpolated is located. Calculate the spatial distance between the pixel to be interpolated and its natural neighbors, and calculate the weight of each natural neighbor based on the distance; the closer the distance, the greater the weight. The linear deformation rate values ​​of each natural neighbor point are weighted and averaged according to their weights to obtain the linear deformation rate interpolation of the pixel to be interpolated. If the Kriging interpolation method is used, the spatial semivariogram model can be calculated and fitted based on the known linear deformation rate of the land surface at the sample points. Using the fitted spatial semi-variogram model, a set of Kriging equations is constructed, and the Kriging weights of each known sample point for the interpolated pixel are obtained by solving the equations. The linear deformation rate values ​​of the land surface at each known sample point are weighted and summed according to the corresponding Kriging weights to obtain the linear deformation rate interpolation and its estimated variance of the pixel to be interpolated. If the inverse distance weighting method is used, a search radius is set with the pixel to be interpolated as the center, and known sample points within the search radius participate in the interpolation; Calculate the distance between the pixel to be interpolated and each sample point participating in the interpolation, and use the reciprocal of the distance raised to the power of p as the weight of the corresponding sample point; The linear deformation rate values ​​of all sample points involved in the interpolation are weighted and averaged according to their weights to obtain the linear deformation rate interpolation value of the pixel to be interpolated.

[0016] By using an iterative optimization algorithm, the unwrapped phase field that minimizes the total weighted cost or the overall "resistance" is found, and the optimized phase unwrapping result is obtained.

[0017] Compared with the prior art, the advantages and positive effects of the present invention are as follows: Leveraging the inherent multivariate scattering characteristics of multipolar SAR data, this study optimizes and corrects the temporal phase unwrapping process for stable scattering point targets. Data constraints are formed using scattering information from different polarization dimensions to correct the inherent defects of traditional single-dimensional phase calculations, systematically eliminating atmospheric delay errors, orbital errors, and residual terrain phase errors. This reduces the impact of various environmental and system interference factors on the temporal phase, controls the accumulation of data bias over long-term observation periods, and improves the computational accuracy of interferometric phase sequences. The phase unwrapping computational logic is standardized, enhancing the authenticity and stability of the scattering phase signal, and providing standardized, high-precision foundational data for subsequent deformation parameter inversion.

[0018] Based on the optimized deformation phase sequence, refined deformation modeling was carried out, and the linear and nonlinear deformation time series of the land surface were separated and extracted to distinguish deformation parameters with different change patterns. Global spatial interpolation processing was performed on deformation data acquired from discrete points to fill spatial gaps between discrete monitoring points. A unified spatial connection standard for deformation data within the region was established to eliminate monitoring blind spots caused by discrete monitoring modes, forming a complete and continuous regional land surface deformation field. The characteristics of local deformation changes and the spatial distribution patterns across the entire region were fully preserved, the spatial expression of land surface deformation was refined, the overall monitoring system for large-scale land surface deformation was improved, and the data output content of deformation monitoring was enriched. Attached Figure Description

[0019] Figure 1 This is a flowchart of the surface deformation monitoring method based on multi-polarization temporal SAR data described in this invention; Figure 2 A flowchart for identifying and extracting stable scattering point targets; Figure 3 A flowchart for optimizing time-series phase unwrapping using multipolar scattering characteristics. Detailed Implementation

[0020] 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 embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0021] In the description of this invention, it should be understood that the terms "length," "width," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientation or positional relationships, are based on the orientation or positional relationships shown in the accompanying drawings and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, in the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0022] See Figure 1 This invention provides a method for monitoring surface deformation based on multi-polarization time-series SAR data, the specific method including: Multi-temporal, multi-polarization synthetic aperture radar (SAR) images of a pre-defined monitoring area under the same orbit were acquired. Radiometric calibration, filtering, and registration preprocessing were performed on all SAR images to obtain a time-registered multi-polarization SAR dataset. Based on this time-registered multi-polarization SAR dataset, stable scatterer point targets within the monitoring area were identified and extracted using the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm. For the extracted stable scatterer point targets, the temporal phase unwrapping of these targets was optimized using multi-polarization scattering characteristics to obtain an optimized deformation phase sequence that removes atmospheric delay errors, orbital errors, and residual topographic phase errors. Geocoding and deformation modeling were performed on the optimized deformation phase sequence to separate and extract the linear and nonlinear deformation time series of the surface in the monitoring area. Spatial interpolation was performed on the extracted linear and nonlinear deformation time series to generate a continuously spatially distributed surface deformation field in the monitoring area, and the deformation monitoring results were output.

[0023] In one embodiment of the present invention, during the data preprocessing stage, all acquired multi-temporal, multi-polarization SAR images are sequentially radiometrically calibrated to convert the original echo signal intensity into backscattering coefficients. Adaptive filtering algorithms are applied to the radiometrically calibrated SAR images of each polarization channel to suppress speckle noise. A SAR image with the best imaging quality in the middle of the time series is used as a common master image, and all other temporal slave images are registered with this common master image. During registration, coarse registration is performed based on orbital parameters and the digital elevation model (DEM), followed by fine registration at the pixel level using cross-correlation or feature matching methods. The image pairs that have completed pixel-level registration are further optimized at the sub-pixel level using spectral diversity to obtain a registration offset field. This registration offset field is then used to resample all slave images, geometrically aligning the SAR images of all temporal phases to form a time-series registered multi-polarization SAR dataset.

[0024] In practical implementation, taking a scenario of monitoring land subsidence in a coastal city area using Sentinel-1 satellite C-band multi-polarization temporal synthetic aperture radar (SAR) imagery as an example, a data preprocessing step is implemented. This involves processing 20 up-orbit, multi-temporal, multi-polarization SAR images covering the monitoring area from the same orbit. These SAR images contain VV and VH dual-polarization data, spanning two years. Radiometric calibration is performed on all SAR images. This process converts the raw echo signal intensity into backscattering coefficients to eliminate sensor system errors. The radiometric calibration process, based on SAR system parameters and calibration constants, converts digital quantization values ​​into physically meaningful values. 0 value.

[0025] In some embodiments, adaptive filtering algorithms are applied to the radiometrically calibrated synthetic aperture radar images of each polarization channel, for example, using a window size of 7. A 7-pixel RefinedLee filter independently filters the VV and VH polarization channels of each image. The filtering process suppresses the inherent speckle noise of the synthetic aperture radar (SAR) images while striving to preserve the texture information of linear features such as building edges and roads. After filtering, the equivalent number of views of the image is improved, which is beneficial for subsequent stable scatterer identification. In practice, one image located in the middle of the time series with the lowest Doppler centroid and lowest atmospheric water vapor content is selected from all 20 temporal SAR images as the common master image. The remaining 19 images are used as slave images and registered with the common master image. The registration process first uses the precise orbital ephemeris data inherent in the SAR images and the publicly available 30-meter resolution digital elevation model (DEM) data for coarse registration to eliminate the main geometric deviations caused by orbital position and terrain undulations. After coarse registration, fine registration is performed at the pixel level using the cross-correlation method.

[0026] In some embodiments, the cross-correlation method uniformly selects hundreds of control points on the common master image and each secondary image, calculates the cross-correlation coefficient of the neighborhood window of each control point, and obtains the pixel-level registration offset by finding the position of the maximum cross-correlation coefficient. When the cross-correlation method fails in some areas due to low coherence, a feature matching method based on SIFT feature points is used to supplement the acquisition of control points. The fine registration process ensures that the registration error between all secondary images and the common master image is less than 0.1 pixels. Optionally, for each pair of synthetic aperture radar images that have completed pixel-level registration, the spectral diversity method is further used for sub-pixel-level registration optimization. The spectral diversity method estimates the sub-pixel-level registration offset of each control point with an accuracy of 0.001 pixels by analyzing the linear phase gradient of the synthetic aperture radar signal in the range and azimuth directions, thereby obtaining a high-precision two-dimensional registration offset field. This offset field finely characterizes the small geometric distortions caused by unmodeled terrain residuals and atmospheric effects between the master and secondary images.

[0027] In practical implementation, the obtained two-dimensional registration offset field is used to resample all 19 secondary images. Resampling employs a sinc function interpolation kernel to precisely align the geometric coordinates of each secondary image to the coordinate grid of the common master image. After resampling, all 20 VV and VH polarimetric synthetic aperture radar (SAR) images from different time phases are geometrically perfectly aligned, forming a time-series strictly registered multi-polarimetric SAR dataset. It can be understood that radiometric calibration normalizes the original signal to comparable backscattering coefficients, adaptive filtering improves the signal-to-noise ratio, and multi-step registration and resampling ensure the geometric consistency of data from different time phases and polarizations, providing a data foundation for constructing a high-precision differential interferometry network. In another implementation scenario, if L-band four-polarimetric SAR data from the ALOS-2 satellite is used, the radiometric calibration parameters and filter window size need to be adjusted accordingly, but the core processing flow remains the same as described above.

[0028] In one embodiment of the present invention, when identifying a stable scattering point target, refer to... Figure 2For a time-registered multi-polarization SAR dataset, the amplitude deviation index and phase stability index of the backscattering coefficient time series under all temporal phases and all polarization channels are calculated for each pixel location. Combining the amplitude deviation index and phase stability index, a multi-dimensional feature space is constructed, and the pixel's position in this multi-dimensional feature space is used as a measure of its scattering stability. Based on this multi-dimensional feature space, a clustering analysis algorithm is used to automatically identify pixels with both low amplitude deviation index and high phase stability index as candidate permanent scatterer points. For the identified candidate permanent scatterer points, their dominant scattering mechanism is distinguished using their scattering matrix information under different polarization channels. Based on a preset scattering mechanism stability threshold, connectivity analysis is performed on candidate permanent scatterer points belonging to the same scattering mechanism. A scattering mechanism stability threshold is set over a time series, and the scattering mechanism time series of each candidate point is calculated and the change in scattering mechanism between adjacent time phases is analyzed. Candidate points whose change is consistently lower than the scattering mechanism stability threshold are identified as permanent scatterer points with stable scattering mechanisms. In the spatial domain, eight-neighbor or four-neighbor connectivity analysis is performed on all permanent scatterer points identified as having stable scattering mechanisms to identify and mark all spatially connected point groups. Isolated point groups or single points with fewer than a preset minimum number of pixels in the point group are considered noise points and are removed. Connected point groups with more than or equal to the preset minimum number of pixels are retained as the final extracted set of stable scatterer point targets.

[0029] When implementing the PS-InSAR algorithm for permanent scatterer synthetic aperture radar interferometry, a differential interferometric network suitable for multi-polarization data is constructed. This network adds interferometric pairs considering polarization basis variations to the conventional short-spatial baseline combination. For each differential interferogram in the network, phase deviations introduced by polarization basis inconsistencies are estimated and corrected using time-registered multi-polarization SAR data on stable scatterer point targets, forming polarization-corrected differential interferometric phases. The terrain phase contribution on stable scatterer point targets is simulated and removed using external digital elevation model (DEM) data to generate residual differential interferometric phases. A graph-based three-dimensional phase unwrapping method is used to jointly unwrap the residual differential interferometric phases in the spatial and temporal domains to obtain the absolute differential interferometric phases. Based on the spatiotemporal correlation of atmospheric phases, a filtering method is used to separate the atmospheric phase components from the absolute differential interferometric phases. The residual digital elevation model error is estimated using spatially continuous and temporally uncorrelated high-frequency phase residuals. By sequentially removing the atmospheric phase component and the residual digital elevation model error from the absolute differential interferometric phase, the optimized residual differential interferometric phase is obtained. This optimized residual differential interferometric phase sequence is the optimized deformation phase sequence.

[0030] In practical implementation, the time-registered multi-polarization synthetic aperture radar dataset covering a coastal city area is processed to identify and extract stable scattering point targets. During the identification process, for each pixel location, the amplitude deviation index and phase stability index of its backscattering coefficient time series under all 20 time phases and the VV and VH polarization channels are calculated. The calculation formula is: ; in: This represents the standard deviation of the backscattering coefficient amplitude of the image point across all time phases and all polarization channels. The mean amplitude is represented by the value of the phase, while the phase stability index is obtained by analyzing the temporal correlation of the phase sequence of all possible differential interferometric pairs for that image point, and characterizes the stability of the phase in the time dimension.

[0031] In some embodiments, a two-dimensional feature space is constructed by combining the calculated amplitude deviation index and phase stability index of each pixel, with the amplitude deviation index as the horizontal axis and the phase stability index as the vertical axis. Each pixel is projected into this feature space, and its position in the feature space, especially its Euclidean distance from the low amplitude deviation and high phase stability region, is used as a comprehensive measure of its scattering stability. Subsequently, the K-means clustering analysis algorithm is used to automatically classify all pixel points in the feature space, with the number of clusters set to 4. After the algorithm converges iteratively, pixels belonging to the "low amplitude deviation index and high phase stability index" combination category are automatically identified as candidate permanent scatterer points.

[0032] In practice, for the tens of thousands of candidate permanent scattering points identified, the dominant scattering mechanism is distinguished by the scattering matrix information under different polarization channels. For each candidate point, the average coherence matrix of all time phases is calculated, and the average coherence matrix is ​​decomposed into eigenvalues. Based on the distribution relationship of the eigenvalues, the dominant scattering mechanism is classified into surface scattering, volume scattering, or secondary scattering. For example, the eigenvalue decomposition results of candidate points extracted from the dihedral structure formed by urban buildings and the ground clearly show that the secondary scattering mechanism is dominant.

[0033] Optionally, connectivity analysis is performed on candidate permanent scatterer points belonging to the same scattering mechanism, based on a preset scattering mechanism stability threshold. The scattering mechanism stability threshold is defined as the scattering mechanism type not being allowed to change between adjacent time phases, and the change in the contribution of the dominant scattering mechanism not exceeding 15%. The scattering mechanism type and contribution time series of each candidate permanent scatterer point across 20 time phases are calculated, and the consistency of the scattering mechanism type and the change in contribution between adjacent time phases are analyzed. Candidate points whose change is consistently lower than the scattering mechanism stability threshold are identified as permanent scatterer points with stable scattering mechanisms. In the spatial domain, eight-neighborhood connectivity analysis is performed on all permanent scatterer points identified as having stable scattering mechanisms.

[0034] It is understandable that connectivity analysis identifies and marks all spatially connected point groups. Isolated point groups or single points with fewer than 5 pixels are considered pseudo-points formed by noise or temporary scatterers and are removed. Connected point groups with 5 or more pixels are retained. These point groups are spatially continuous and usually correspond to real and reliable stable targets such as building roofs, bridges, and stable radar corner reflectors. All pixels in these connected point groups are used as the final set of stable scatterer point targets for subsequent differential interferometry processing. In another implementation scenario of landslide monitoring in mountainous areas, the dominant scattering mechanism may be surface scattering, and the neighborhood window and minimum point group size threshold of connectivity analysis can be adjusted according to the actual terrain fragmentation, but the core discrimination logic is consistent with the above description.

[0035] In practice, the PS-InSAR algorithm for permanent scatterer synthetic aperture radar interferometry is implemented based on the final extracted set of stable scatterer point targets. A differential interferometric network suitable for multi-polarization data is constructed. The differential interferometric network adds interferometric pairs that consider polarization basis changes, based on the conventional combination of short spatial baselines (less than 150 meters) and short temporal baselines (less than 60 days). For example, data from the same time phase but different polarization channels (VV and VH) are combined to form polarization interferometric pairs. The entire differential interferometric network consists of 45 differential interferograms. For each differential interferogram in the differential interferometric network, the phase deviation introduced by the inconsistency of polarization basis is estimated and corrected on the stable scatterer point targets using the multi-polarization information of the multi-polarization synthetic aperture radar data with time-registration. For each stable scatterer point, a polarization phase difference is estimated and compensated by comparing its scattering vectors under different polarization channels in the master and slave images, forming a polarization-corrected differential interferometric phase.

[0036] In some embodiments, the terrain phase contribution on stable scatterer point targets is simulated and removed using external 30-meter resolution digital elevation model (DEM) data to generate residual differential interferometric phase. A graph-based three-dimensional phase unwrapping method is used on the residual differential interferometric phase to jointly unwrap the phase in the spatial and temporal domains to obtain the absolute differential interferometric phase. Based on the spatiotemporal correlation of atmospheric phase, the atmospheric phase component is separated from the absolute differential interferometric phase using a combination of high-pass filtering (in the time domain) and low-pass filtering (in the spatial domain). The residual digital elevation model error is estimated using spatially continuous and temporally uncorrelated high-frequency phase residuals. The atmospheric phase component and the residual digital elevation model error are sequentially removed from the absolute differential interferometric phase to obtain the optimized residual differential interferometric phase. The sequence consisting of the optimized residual differential interferometric phases of all 45 interferograms is the optimized deformation phase sequence.

[0037] In one embodiment of the present invention, when optimizing time-series phase unwrapping using multipolar scattering characteristics, see [reference needed]. Figure 3 For each stable scattering point target, its scattering matrix or coherence matrix is ​​extracted across all time phases and polarization channels. Eigenvalue decomposition is performed on the scattering matrix or coherence matrix to obtain the dominant scattering mechanism and its corresponding scattering contribution of the stable scattering point target in different time phases. Based on the time-series stability of the dominant scattering mechanism and its corresponding scattering contribution, the reliability of the stable scattering point target's phase in different interferometric pairs is evaluated, and a polarization-based weighting coefficient is assigned to the phase of each interferometric pair. During phase unwrapping, the polarization-based weighting coefficient is introduced into the unwrapping smoothing constraint or cost function, defining a phase unwrapping cost function. This cost function contains the weighted sum of squares of the differences between the unwrapped phase gradient and the wrapped phase gradient. The polarization-based weighting coefficient is used as the weighting factor for the corresponding interferometric pair and the corresponding point target in the cost function. When using the cost function to minimize phase unwrapping, the phase gradient difference of the high-weighted interferometric pair contributes more to the total cost. An iterative optimization algorithm is used to find the unwrapped phase field that minimizes the weighted total cost, making the phase unwrapping process more inclined to follow the change path of high-weighted, high-reliability phases. For the unwrapped absolute phase, additional constraints are constructed using the phase relationship between the multi-polarization channels to perform cross-validation and consistency correction on the unwrapping results, thereby optimizing the phase unwrapping results.

[0038] In practice, based on a set of stable scattering point targets covering coastal urban areas, the scattering matrix of each stable scattering point target is extracted across all 20 time phases and the two polarization channels, VV and VH. For each time phase, the scattering matrix of each stable scattering point target is a 2x2 complex matrix containing co-polarization and cross-polarization information. The coherence matrix of each stable scattering point target in each time phase is calculated. By performing eigenvalue decomposition on the coherence matrix, the dominant scattering mechanism and corresponding scattering contribution of the stable scattering point target in different time phases are obtained. Eigenvalue decomposition produces three eigenvalues ​​and their corresponding eigenvectors. The eigenvector corresponding to the largest eigenvalue indicates the dominant scattering mechanism, while the scattering contribution is represented by the percentage of each eigenvalue to the sum, thus forming a time series of scattering mechanism and contribution of each point target across all time phases.

[0039] In some embodiments, the reliability of the phase of a stable scattering point target in different interferometric pairs is evaluated based on the time-series stability of the dominant scattering mechanism and the corresponding scattering contribution. The consistency of the dominant scattering mechanism and the magnitude of the scattering contribution variation for each stable scattering point target across the two time phases involved in each interferogram pair are calculated. Interferometric pairs with smaller scattering contribution variations and consistent scattering mechanisms are considered to have more reliable phases. Based on this evaluation, a polarization-based weighting coefficient is assigned to each stable scattering point target in each interferometric pair. The weighting coefficient ranges from 0 to 1, with higher values ​​indicating higher phase reliability. Table 1 shows the evaluation results of the polarization weighting coefficients for three stable scattering point targets in some interferometric pairs. Table 1: Polarization weighting coefficients of stable scattering point targets in some interferometric pairs

[0040] In practical implementation, during the phase unwrapping process, weighting coefficients based on polarization information are introduced into the cost function of unwrapping. The cost function of phase unwrapping is defined as the sum of squared weighted gradient differences of all stable scattering point targets over all interference pairs and spatial neighborhoods. The expression is: ; in: This represents the total number of stable scattering point targets. This represents the total number of differential interferograms. Let represent the set of point targets that are spatially adjacent to the i-th point target. These are the weighting coefficients based on polarization information assigned to the i-th point target in the j-th differential interferogram. This represents the gradient of the unwrapped phase in the spatial neighborhood. The gradient of the entangled phase in the same spatial neighborhood is represented by the weight coefficients based on polarization information, which serve as weight factors for the corresponding interference pairs and point targets in the cost function. The weight factors for high-reliability phases are large, while those for low-reliability phases are small. When using the minimization cost function for phase unwrapping, the phase gradient difference of high-weighted interference pairs contributes more to the total cost, forcing the unwrapping result to be more inclined to satisfy the continuity of high-reliability phases. Through iterative optimization algorithms such as weighted least squares or network flow algorithms, the unwrapped phase field that minimizes the weighted total cost is found, and the optimized phase unwrapping result is obtained.

[0041] Optionally, for the unwrapped absolute phase, additional constraints are constructed using the phase relationship between the multi-polarization channels to perform cross-validation and consistency correction on the unwrapping results. For each stable scattering point target, its unwrapped absolute phase in the VV and VH polarization channels should theoretically satisfy a specific phase relationship under noise-free and ideal scattering conditions. For example, on a point target dominated by secondary scattering, the dual-polarization phase difference should be close to 0 or π. Based on this physical relationship, a consistency check threshold is set. For point targets that do not satisfy the threshold relationship, the unwrapped phase path in its local region is checked, and the integer ambiguity of its unwrapped phase is adjusted iteratively to restore the consistency of the phase relationship between the multi-polarization channels, thereby correcting unwrapping errors and reducing the uncertainty of phase unwrapping.

[0042] It is understandable that by introducing weighting coefficients based on polarization scattering stability, the phase unwrapping process can distinguish the reliability of phase information for different point targets and different interferometric pairs, and give higher priority to high-reliability phases in the unwrapping process. This helps to obtain more reliable unwrapping results in areas with drastic phase gradient changes or high noise. In another implementation scenario, if fully polarimetric synthetic aperture radar data is used, richer polarization information can be used to construct more refined weighting coefficients. However, the core method of incorporating weighting coefficients into the phase unwrapping cost function is consistent with the above description.

[0043] In one embodiment of the present invention, when geocoding and modeling the optimized deformation phase sequence, the orbital parameters and imaging geometry of the radar system are used to convert the optimized deformation phase sequence from the slant range geometry of the radar coordinate system to the ground geometry in the geographic coordinate system, thus completing the geocoding. Assuming that the deformation consists of a linear trend term, a seasonal periodic term, and an irregular nonlinear residual term, a time series model of surface deformation is constructed. Using the geocoded optimized deformation phase sequence as observations, the parameters of the time series model of surface deformation are solved using least squares estimation or singular value decomposition. The coefficients of the linear trend term are directly extracted from the solved model parameters as the surface linear deformation rate. Subtracting the fitted linear trend term and the known seasonal periodic term from the optimized deformation phase sequence yields the nonlinear deformation time series.

[0044] In practice, the optimized deformation phase sequence is geocoded and modeled. Using precise orbital ephemeris data and imaging geometry parameters from Sentinel-1 satellite imagery, including satellite position, velocity vector, slant range, and radar wavelength, the optimized deformation phase sequence of each stable scatterer point target is transformed from its radar slant range coordinate system to latitude, longitude, and elevation coordinates in the WGS84 geographic coordinate system. This transformation process is achieved through rigorous calculations using the radar range-Doppler equation and the Earth ellipsoid model. After geocoding, the deformation phase value of each point target is associated with its precise geodetic coordinates.

[0045] In some embodiments, it is assumed that surface deformation over time consists of a linear trend term, a seasonal periodic term, and an irregular nonlinear residual term. A time series model of surface deformation is constructed, where the linear trend term represents the long-term average rate of deformation, the seasonal periodic term characterizes periodic deformation caused by factors such as groundwater extraction and thermal expansion, and the nonlinear residual term includes unmodeled deformation signals, noise, and possible abrupt deformations. For a point target, the time series model is constructed. Deformation phase Its time series model is expressed as: ; in: It is a linear deformation rate coefficient, which is directly related to the annual average deformation rate in the radar line-of-sight direction. and It is the amplitude coefficient of the seasonal periodic term, and f is the frequency of the seasonal variation, usually taken as 1 year. -1 , It is a nonlinear residual term. Let represent the annualized time of the acquisition time of the k-th image relative to the start time of the time series. Using the geocoded optimized deformation phase sequence as the observations, the parameters of the above surface deformation time series model, i.e., the coefficient vector, are solved using the least squares estimation method. .

[0046] In specific implementation, taking five representative stable scatterer point targets in the monitoring area as an example, the least squares estimation results of their deformation model parameters are shown in the table below. From the solved model parameters, the coefficient α of the linear trend term is directly extracted as the surface linear deformation rate of each point target. Positive values ​​indicate deformation away from the satellite, and negative values ​​indicate deformation towards the satellite. See Table 2.

[0047] Table 2: Parameter Estimation Results of Stable Scattering Point Target Deformation Model

[0048] Optionally, when using least squares estimation to solve for model parameters, for large-scale calculations involving tens of thousands of point targets across the entire monitoring area, block matrix operations or iterative optimization algorithms can be employed to improve computational efficiency. In another implementation scenario, if the deformation time series exhibits a significant nonlinear trend or complex periodic components, singular value decomposition (SVD) can be used instead of least squares estimation to solve for the parameters of the surface deformation time series model. SVD can more stably handle ill-conditioned design matrices and obtain the minimum norm solution for the model parameters. The model parameters are then subtracted from the optimized deformation phase sequence. , , The linear trend term and seasonal periodic term fitted to the known frequency f, and the remaining part constitute the nonlinear deformation time series of each point target. Nonlinear deformation time series contain unmodeled deformation signals, residual atmospheric and orbital noise, and information on possible abrupt deformations, such as short-period deformations caused by local geological activity or engineering construction.

[0049] It is understandable that geocoding maps deformation information from radar geometric space to real-world geographic space, which is the basis for the visualization and analysis of deformation results. Time series modeling and parameter solving decompose the mixed phase observations into linear deformation rates, periodic deformations, and nonlinear residuals with clear physical meanings. The linear deformation rate directly reflects the long-term trend of surface deformation and is the core result of subsidence monitoring. The nonlinear deformation time series provides a data basis for further analysis of transient or sporadic deformation events. In another implementation scenario, if there is no obvious seasonal deformation in the monitoring area, the time series model can contain only linear terms and nonlinear residual terms. Its parameter solving and separation process is essentially the same as described above.

[0050] In one embodiment of the present invention, when spatially interpolating the extracted linear deformation rate and nonlinear deformation time series of the land surface, the geographical location of the extracted stable scatterer point target is used as a known sample point, and its corresponding linear deformation rate is used as a sample value. Based on the natural neighbor method, Kriging interpolation, or inverse distance weighting, spatial interpolation calculation is performed on the linear deformation rate of the known sample points. If the natural neighbor method is used, for each pixel to be interpolated within the monitoring area, its natural neighbors are found among the known sample points. The spatial distance between the pixel to be interpolated and each of its natural neighbors is calculated, and the weight of each natural neighbor is calculated based on the distance. The linear deformation rate values ​​of each natural neighbor are then weighted and averaged according to their weights to obtain the linear deformation rate interpolation value of the pixel to be interpolated. If the Kriging interpolation method is used, the spatial semivariogram model is calculated and fitted based on the linear deformation rate values ​​of the known sample points, and the fitted spatial semivariogram is used... A semi-variogram model is used to construct and solve the Kriging equations to obtain the Kriging weights for each known sample point to be interpolated. The linear deformation rate values ​​of each known sample point are then weighted and summed according to their corresponding Kriging weights to obtain the linear deformation rate interpolation and its estimated variance for the pixel to be interpolated. Alternatively, an inverse distance weighting method is used, with a search radius centered on the pixel to be interpolated. The distance between the pixel to be interpolated and each sample point within the search radius is calculated, and the reciprocal of the p-th power of the distance is used as the weight. The linear deformation rate values ​​of the sample points are then weighted and averaged according to their weights to obtain the interpolation, generating a continuous spatial distribution raster map of linear deformation rate covering the entire monitoring area. For each time point in the nonlinear deformation time series, the nonlinear deformation value of the stable scattering point target at the corresponding time point is used as the sample value. Using the same spatial interpolation method as for linear deformation rate, independent spatial interpolation is performed on the nonlinear deformation sample at each time point, generating a nonlinear deformation spatial distribution raster map corresponding to each time point. By integrating the continuous spatial distribution of linear deformation rate raster map with the nonlinear deformation spatial distribution raster map at all time points, a complete surface deformation field is formed.

[0051] In practice, based on the geographical location of stable scattering point targets within the monitoring area and their corresponding linear deformation rates of the land surface, spatial interpolation is performed to generate a continuously spatially distributed land surface deformation field. The total number of stable scattering point targets extracted within the monitoring area exceeds 50,000. The geographical coordinates of these point targets and their corresponding linear deformation rate values ​​constitute a known sample point set. The interpolation objective is to assign a linear deformation rate value to each pixel in the regular geographic grid covering the entire monitoring area. The spatial resolution of the regular geographic grid is set to 50 meters.

[0052] In some embodiments, spatial interpolation of the surface linear deformation rate of known sample points is performed based on the natural neighbor method. For each pixel to be interpolated within the monitoring area, its natural neighbors are found among the known sample points. This search is achieved by constructing a Deloni triangulation network of all known sample points. A natural neighbor is a sample point in the Deloni triangulation network that is directly adjacent to the Voronoi polygon cell containing the pixel to be interpolated. The spatial distance between the pixel to be interpolated and each of its natural neighbors is calculated, and the weight of each natural neighbor is calculated based on the distance. The weight calculation adopts the Sibyl interpolation weight of the natural neighbor method, which is the proportion of the area of ​​the pixel to be interpolated falling into the Voronoi cell defined by the natural neighbors to the area falling into the original Voronoi cell defined by all sample points. The surface linear deformation rate value of each natural neighbor is weighted and averaged according to its corresponding Sibyl weight to obtain the linear deformation rate interpolation result of the pixel to be interpolated. The calculation process is expressed as the formula: ; in: This represents the linear deformation rate interpolation result of the pixel to be interpolated, where n represents the number of natural neighbors. The Sibyl interpolation weight of the i-th natural neighbor is determined by the area ratio and satisfies the following conditions: , This represents the linear deformation rate value of the land surface corresponding to the i-th natural neighbor. Iterate through all the pixels to be interpolated within the monitoring area and repeat the above process to generate a linear deformation rate raster map covering the entire monitoring area with continuous spatial distribution.

[0053] In practice, for each time point of the nonlinear deformation time series, the nonlinear deformation value of the stable scattering point target at the corresponding time point is used as the sample value. The nonlinear deformation time series contains twenty time points, corresponding to the acquisition times of twenty synthetic aperture radar images. The same natural neighbor spatial interpolation method as the linear deformation rate is used to perform independent spatial interpolation on the nonlinear deformation sample at each time point. For each time point, based on the nonlinear deformation values ​​of all stable scattering point targets at that time, the Sibyl interpolation weight of each pixel to be interpolated relative to its current natural neighbor is recalculated and a weighted average is performed to generate a nonlinear deformation spatial distribution raster map corresponding to each time point. Finally, twenty nonlinear deformation spatial distribution raster maps are obtained.

[0054] Optionally, in another implementation scenario, if the known sample points are spatially unevenly distributed and interpolation uncertainty needs to be evaluated, the Kriging interpolation method is used. Based on the known surface linear deformation rate values ​​of the sample points, a spatial semi-variogram model is calculated and fitted. The semi-variogram model can be an exponential model or a Gaussian model. The fitted spatial semi-variogram model is used to construct and solve a system of Kriging equations to obtain the Kriging weights for each known sample point to be interpolated. The surface linear deformation rate values ​​of each known sample point are then weighted and summed according to their corresponding Kriging weights to obtain the interpolated pixel. In another implementation scenario, if computational resources are limited and high interpolation efficiency is required, the inverse distance weighting method is used for the linear deformation rate interpolation of pixels and its estimated variance. A search radius of 500 meters is set with the pixel to be interpolated as the center. Known sample points within the search radius participate in the interpolation. The distance between the pixel to be interpolated and each sample point participating in the interpolation is calculated. The reciprocal of the distance raised to the power of p is used as the weight of the corresponding sample point. The surface linear deformation rate values ​​of all sample points participating in the interpolation are weighted and averaged according to their weights to obtain the linear deformation rate interpolation of the pixel to be interpolated.

[0055] It is understandable that integrating the linear deformation rate raster map with the nonlinear deformation spatial distribution raster map at all time points forms a complete surface deformation field. The linear deformation rate raster map characterizes the long-term spatial trend of deformation, while the nonlinear deformation raster map sequence reveals the dynamic details of deformation in the time dimension. The integrated surface deformation field is continuous in space and contains both linear trends and nonlinear details in time, which can comprehensively reflect the surface deformation characteristics of the monitoring area.

[0056] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for monitoring surface deformation based on multi-polarization temporal SAR data, characterized in that, Includes the following steps: Collect multi-temporal and multi-polarization synthetic aperture radar (SAR) images of the preset monitoring area and the same orbit, and perform radiometric calibration, filtering and registration preprocessing on all SAR images to obtain a time-registered multi-polarization SAR dataset. Based on the time-registered multipolar SAR dataset, stable scatterer point targets within the monitoring area are identified and extracted using the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm. For the extracted stable scatterer point targets, the temporal phase unwrapping of the stable scatterer point targets is optimized by utilizing the multi-polarization scattering characteristics to obtain an optimized deformation phase sequence that removes atmospheric delay error, orbital error and residual terrain phase error; Geocoding and deformation modeling are performed on the optimized deformation phase sequence to separate and extract the surface linear deformation rate and nonlinear deformation time series of the monitoring area; Spatial interpolation is performed on the extracted linear and nonlinear deformation time series of the land surface to generate a continuously spatially distributed land surface deformation field in the monitoring area, and the deformation monitoring results are output. The process of performing radiometric calibration, filtering, and registration preprocessing on all SAR images to obtain a time-registered multi-polarization SAR dataset is as follows: All acquired multi-temporal, multi-polarization SAR images are sequentially radiometrically calibrated to convert the original echo signal intensity into backscattering coefficients, thereby eliminating sensor system errors. For the radiometrically calibrated SAR images of each polarization channel, an adaptive filtering algorithm is applied to suppress the speckle noise of the images while preserving the edge and texture information of the images. A common master image is selected as the SAR image with the best imaging quality in the middle of the time series, and all other secondary images in the time series are registered with the common master image. During registration, coarse registration is performed based on orbital parameters and the digital elevation model (DEM), and fine registration is performed at the pixel level using cross-correlation or feature matching methods to ensure that the registration error between all slave images and the common master image is less than a preset threshold. For image pairs that have completed pixel-level registration, the spectral diversity method is further used to optimize the sub-pixel-level registration and obtain the registration offset field. All images are resampled using the registration offset field to geometrically align SAR images of all time phases, forming a time-registered multi-polarization SAR dataset.

2. The surface deformation monitoring method based on multi-polarization temporal SAR data according to claim 1, characterized in that, Based on the time-registered multi-polarization SAR dataset, stable scatterer point targets within the monitoring area are identified and extracted using the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm. Specifically: For the time-registered multi-polarization SAR dataset, for each pixel location, the amplitude deviation index and phase stability index of its backscattering coefficient time series under all time phases and all polarization channels are calculated; Combining the amplitude deviation index and the phase stability index, a multidimensional feature space is constructed, and the position of the pixel in the multidimensional feature space is used as a measure of its scattering stability. Based on the multidimensional feature space, a clustering analysis algorithm is used to automatically identify pixels that simultaneously have low amplitude deviation index and high phase stability index as candidate permanent scattering points. For the identified candidate permanent scatterers, their dominant scattering mechanisms are distinguished by their scattering matrix information under different polarization channels. By combining the preset stability threshold of scattering mechanism, connectivity analysis is performed on candidate permanent scatterer points belonging to the same scattering mechanism to eliminate isolated noise points, and finally extract the set of stable scatterer point targets that are spatially coherently distributed within the monitoring area.

3. The surface deformation monitoring method based on multi-polarization temporal SAR data according to claim 2, characterized in that, The permanent scatterer synthetic aperture radar interferometry (PS-InSAR) algorithm specifically includes: A differential interferometric network suitable for multipolar data is constructed. The differential interferometric network adds interferometric pairs that consider polarization basis changes on the basis of conventional short-spatial-temporal baseline combinations. For each differential interferogram in the differential interferometry network, on the stable scatterer point target, the phase deviation introduced by the polarization basis inconsistency is estimated and corrected using the multi-polarization information in the time-registered multi-polarization SAR dataset to form the polarization-corrected differential interferometric phase. The terrain phase contribution on the stable scatterer point target is simulated and removed using external digital elevation model (DEM) data to generate residual differential interferometric phase. For the residual differential interference phase, a graph-based three-dimensional phase unwrapping method is used to jointly unwrap the phase in the spatial and temporal domains to obtain the absolute differential interference phase; Based on the spatiotemporal correlation of atmospheric phase, atmospheric phase components are separated from the absolute differential interferometric phase using a filtering method. The residual digital elevation model error is estimated by using spatially continuous and temporally uncorrelated high-frequency phase residuals. The atmospheric phase component and the residual digital elevation model error are sequentially removed from the absolute differential interferometric phase to obtain the optimized residual differential interferometric phase. The optimized residual differential interferometric phase sequence is the optimized deformation phase sequence.

4. The method for monitoring surface deformation based on multi-polarization temporal SAR data according to claim 1, characterized in that, The step of optimizing the temporal phase unwrapping of the extracted stable scattering point target using multi-polarization scattering characteristics is as follows: For each stable scattering point target, extract its scattering matrix or coherence matrix under all time phases and all polarization channels; The scattering matrix or coherence matrix is ​​decomposed into eigenvalues ​​to obtain the dominant scattering mechanism and corresponding scattering contribution of the stable scattering point target in different time phases. Based on the dominant scattering mechanism and the time-series stability of the corresponding scattering contribution, the reliability of the phase of the stable scattering point target in different interferometric pairs is evaluated, and a weighting coefficient based on polarization information is assigned to the phase of each interferometric pair. During the phase unwrapping process, the weighting coefficients based on polarization information are introduced into the smoothing constraint or cost function of unwrapping, so that the phase unwrapping process tends to follow the change path of high-weight, high-reliability phases. For the unwrapped absolute phase, additional constraints are constructed using the phase relationship between the multi-polarization channels. Cross-validation and consistency correction are performed on the unwrapping results to correct unwrapping errors, thereby optimizing the phase unwrapping results and reducing the uncertainty of phase unwrapping.

5. The method for monitoring surface deformation based on multi-polarization temporal SAR data according to claim 1, characterized in that, The optimized deformation phase sequence is geocoded and modeled for deformation, and the linear deformation rate and nonlinear deformation time series of the monitored area are separated and extracted, specifically as follows: Using the orbital parameters and imaging geometry of the radar system, the optimized deformation phase sequence is converted from the slant range geometry of the radar coordinate system to the ground geometry in the geographic coordinate system, thus completing the geocoding. Assuming that deformation consists of a linear trend term, a seasonal periodic term, and an irregular nonlinear residual term, a time series model of surface deformation is constructed. Using the geocoded optimized deformation phase sequence as observations, the parameters of the time series model of the surface deformation are solved by least squares estimation or singular value decomposition. The coefficients of the linear trend term are directly extracted from the solved model parameters and used as the linear deformation rate of the land surface. Subtracting the fitted linear trend term and the known seasonal periodic term from the optimized deformation phase sequence, the remaining part is the nonlinear deformation time series, which includes unmodeled deformation signals, residual noise, and abrupt deformation information.

6. The method for monitoring surface deformation based on multi-polarization temporal SAR data according to claim 1, characterized in that, The step of spatially interpolating the extracted linear and nonlinear deformation time series of the land surface to generate a continuously spatially distributed land surface deformation field in the monitoring area is as follows: The geographical location of the extracted stable scatterer point target is used as the known sample point, and the corresponding surface linear deformation rate is used as the sample value. Based on the natural neighbor method, Kriging interpolation method or inverse distance weighting method, spatial interpolation calculation is performed on the surface linear deformation rate of the known sample points to generate a continuous spatial distribution of linear deformation rate raster map covering the entire monitoring area. For each time point in the nonlinear deformation time series, the nonlinear deformation value of the stable scattering point target at the corresponding time point is used as a sample value; Using the same spatial interpolation method as the linear deformation rate, independent spatial interpolation is performed on the nonlinear deformation sample at each time point to generate a nonlinear deformation spatial distribution raster map corresponding to each time point. The linear deformation rate raster map with continuous spatial distribution is integrated with the nonlinear deformation spatial distribution raster map at all time points to form a complete surface deformation field. The surface deformation field is continuous in space and contains linear trends and nonlinear details in time.

7. The method for monitoring surface deformation based on multi-polarization temporal SAR data according to claim 2, characterized in that, The process involves combining a preset scattering mechanism stability threshold with connectivity analysis of candidate permanent scatterer points belonging to the same scattering mechanism to eliminate isolated noise points. Specifically: A stability threshold for the scattering mechanism is set over a time series, which defines the maximum allowable range of variation of the scattering mechanism between adjacent time phases; Calculate the time series of scattering mechanisms for each candidate permanent scatterer point and analyze the changes in scattering mechanisms between adjacent time phases; Candidate points whose change in magnitude is consistently below the stability threshold of the scattering mechanism are identified as permanent scattering points with a stable scattering mechanism. In the spatial domain, connectivity analysis of eight-neighbor or four-neighbor domains is performed on all permanent scattering points determined to have stable scattering mechanisms. Identify and mark all spatially connected point groups, and discard isolated point groups or single points with fewer than a preset minimum number of pixels as noise points. The group of connected points with a number of pixels greater than or equal to a preset lower limit is retained as the final extracted stable scattering point target.

8. The method for monitoring surface deformation based on multi-polarization temporal SAR data according to claim 4, characterized in that, During phase unwrapping, the weighting coefficients based on polarization information are introduced into the smoothing constraint or cost function of the unwrapping process, specifically as follows: Define a cost function for phase unwrapping, which includes the weighted sum of squares of the differences between the unwrapped phase gradient and the wrapped phase gradient; The weight coefficients based on polarization information are used as weight factors for the corresponding interference pairs and corresponding point targets in the cost function. The weight factors for high reliability phases are large, and the weight factors for low reliability phases are small. When performing phase unwrapping using the cost function minimization, high-weighted interference contributes more to the phase gradient difference in the total cost, forcing the unwrapping result to be more inclined to satisfy the continuity of the high-reliability phase. By using an iterative optimization algorithm, the unwrapped phase field that minimizes the total weighted cost is found, and the optimized phase unwrapping result is obtained.

9. The method for monitoring surface deformation based on multi-polarization temporal SAR data according to claim 6, characterized in that, The spatial interpolation calculation of the surface linear deformation rate of the known sample points based on the natural neighbor method, Kriging interpolation method, or inverse distance weighting method is specifically as follows: If the natural neighbor method is used, for each pixel to be interpolated in the monitoring area, its natural neighbor is found among the known sample points. The natural neighbor is the sample point in the Thiessen polygon that is adjacent to the polygon in which the pixel to be interpolated is located. Calculate the spatial distance between the pixel to be interpolated and its natural neighbors, and calculate the weight of each natural neighbor based on the distance; the closer the distance, the greater the weight. The linear deformation rate values ​​of each natural neighbor point are weighted and averaged according to their weights to obtain the linear deformation rate interpolation of the pixel to be interpolated. If the Kriging interpolation method is used, the spatial semivariogram model can be calculated and fitted based on the known linear deformation rate of the land surface at the sample points. Using the fitted spatial semi-variogram model, a set of Kriging equations is constructed, and the Kriging weights of each known sample point for the interpolated pixel are obtained by solving the equations. The linear deformation rate values ​​of the land surface at each known sample point are weighted and summed according to the corresponding Kriging weights to obtain the linear deformation rate interpolation and its estimated variance of the pixel to be interpolated. If the inverse distance weighting method is used, a search radius is set with the pixel to be interpolated as the center, and known sample points within the search radius participate in the interpolation; Calculate the distance between the pixel to be interpolated and each sample point participating in the interpolation, and use the reciprocal of the distance raised to the power of p as the weight of the corresponding sample point; The linear deformation rate values ​​of all sample points involved in the interpolation are weighted and averaged according to their weights to obtain the linear deformation rate interpolation of the pixel to be interpolated. By using an iterative optimization algorithm, the unwrapped phase field that minimizes the total weighted cost or the overall "resistance" is found, and the optimized phase unwrapping result is obtained.

Images