A method and system for atmospheric delay correction of interferometric synthetic aperture radar

Abstract

This invention discloses an atmospheric delay correction method and system for interferometric synthetic aperture radar (InSAR), belonging to the field of radar signal correction technology. It includes: performing time-series analysis on a preprocessed SAR image set to obtain multiple interferograms; performing phase unwrapping to obtain multiple unwrapped phases; decomposing the observation matrix composed of the multiple unwrapped phases to obtain multiple independent components; removing independent components related to vertical layer delay from the observation matrix based on the correlation coefficient and periodicity coefficient between each independent component and terrain elevation to obtain a residual phase matrix; and performing turbulent delay correction on the residual phase matrix to obtain a corrected interferogram phase set deformation time series. This invention effectively suppresses the interference of atmospheric delay errors on deformation signals, preserves the true deformation information to the greatest extent, and significantly improves the accuracy of InSAR deformation monitoring.

Description

Technical Field

[0001] This invention relates to the field of radar signal correction technology, and in particular to an atmospheric delay correction method and system for interferometric synthetic aperture radar. Background Technology

[0002] In recent years, Interferometric Synthetic Aperture Radar (InSAR) has become an effective technique for detecting and mapping ground deformation caused by natural activities or human influences. It can monitor deformation over large areas with accuracy ranging from centimeters to millimeters. However, atmospheric variations affect the propagation of radar signals as they pass through the troposphere during two separate SAR acquisitions. Factors such as temperature, pressure, and humidity change over time and space, altering the speed and propagation of radar signals. This results in atmospheric delay, significantly impacting the accuracy of InSAR measurements. This delay can introduce errors in deformation estimation, potentially reaching tens of centimeters. Therefore, the impact of atmospheric delay cannot be ignored to ensure the accuracy of InSAR measurements.

[0003] The existing technology mainly adopts the correction method based on interferometric phase. The implementation process is as follows: by constructing an empirical relationship between the interferometric phase and the terrain elevation (such as a linear model or a power law model), the model parameters are obtained by least squares fitting, and then the vertical stratification delay component estimated by the model is subtracted from the interferometric phase to achieve atmospheric delay correction.

[0004] However, the correction method based on interferometric phase uses a fixed form of parameterized model (linear or power law), which cannot adapt to the complex nonlinear relationship between topography and vertical stratification delay in areas with significant elevation changes. This results in insufficient vertical stratification delay correction, a large impact from residual atmospheric errors, and reduced accuracy of InSAR deformation monitoring. Summary of the Invention

[0005] Therefore, it is necessary to provide an atmospheric delay correction method and system for interferometric synthetic aperture radar to address the aforementioned technical problems.

[0006] This invention provides an atmospheric delay correction method for interferometric synthetic aperture radar, comprising: Acquire Sentinel-1 SAR data captured by interferometric synthetic aperture radar and arrange them according to the capture time sequence to obtain a preprocessed SAR image set; StaMPS was used to perform time series analysis on the preprocessed SAR image set to obtain multiple interferograms; and phase unwrapping was performed on the multiple interferograms to obtain multiple unwrapped phases; The observation matrix is ​​composed of multiple unwrapped phases. The observation matrix is ​​decomposed by variational Bayesian independent component analysis to distinguish the components representing the layered delay from the independent components by utilizing the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components. The independent components related to vertical stratification delay are identified based on the correlation coefficient and periodicity coefficient between each independent component and the terrain elevation. The independent components related to vertical stratification delay are then removed from the observation matrix to obtain the residual phase matrix. Turbulent delay correction was performed on the residual phase matrix using phase stacking and spatiotemporal filtering methods to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar.

[0007] Optionally, phase unwrapping is performed on multiple interferograms to obtain multiple unwrapped phases, specifically including: Each interferogram is considered as a linear combination of the propagating phase differences between various temporal and spatial variation sources, based on the following formula: ; in, To untangle the phase, For the deformed phase component, For orbital error phase, For residual topographic phase, It is random noise. For the scattering phase, This refers to the tropospheric phase.

[0008] Optionally, the observation matrix is ​​composed of multiple unwrapped phases, and then decomposed using variational Bayesian independent component analysis to distinguish between components representing the layered delay and independent components based on the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components, specifically including: The observation matrix is ​​decomposed using variational Bayesian independent component analysis based on the following formula: ; in, For the observation matrix, A It is a mixed matrix. S For the reconstructed source matrix, N Assuming it is zero-mean Gaussian noise; The mixing matrix and the reconstructed source matrix are described as the vbICA model based on the following formula: ={ A , Λ , S , q , θ}; θ ={ Π ,μ , β}; in, A collection of variables and parameters. Λ For Gaussian noise accuracy, q For discrete variables, μ The center location or average value of each Gaussian component. β The reciprocal of the variance of each Gaussian component. Π The weight or probability of each Gaussian component; Given the observation matrix, calculate Posterior probability: ; ; ; ; ; in, This represents the true posterior probability. As a prior distribution, As the normalization factor, For the variational approximation posterior of each parameter, for The first in One variable, For joint probability, for X The marginal likelihood logarithm, F negative free energy ， for Entropy; Multiple independent components are constructed based on the mixing matrix and the reconstructed source matrix.

[0009] Optionally, the correlation coefficient between each independent component and topographic elevation can be determined based on the following formula: ; in, For the first k The correlation coefficient of each independent component For the first k The first of the independent components i Spatial relationship of pixels, For the first k The mean of the spatial relationships of all pixels of each independent component. H i For the first i The elevation estimated by each InSAR pixel This is the average elevation of all pixels; The periodicity coefficients between each independent component and the terrain elevation are determined based on the following formula: ; in, For the first k Power spectrum of time series of individual components F For Fourier transform, For the first k The time series of 3 independent components, where ()* represents complex conjugation. L This represents the length of the time series vector.

[0010] Optionally, the residual phase matrix is ​​corrected for turbulence delay using a phase stacking method and a spatiotemporal filtering method to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar, specifically including: The atmospheric delay for each acquisition date in the residual phase matrix is ​​estimated using the phase stacking method based on the following equation: ; The atmospheric delay for each SAR date is estimated iteratively based on the atmospheric noise figure using the following formula: ; = ; in, and This refers to the time sequence number of the SAR image. N The number of SAR image pairs with the same time interval. For the first i Using the background as a reference, the phase difference between consecutive image pairs constructed under the same time interval conditions. For the first j Using the background as a reference, the phase difference between consecutive image pairs constructed under the same time interval conditions. For the first i Scene Pixel m Atmospheric phase delay at that location The average atmospheric phase for all pixels. M This represents the total number of pixels. The RMS value of the atmospheric phase with the highest noise level; The interferograms are sorted according to the magnitude of the atmospheric phase in each SAR acquisition. Starting from the date with the highest noise level, the atmospheric delay before that date is used to correct the interferograms after that date, resulting in the phase and deformation time series corrected by the interferometric synthetic aperture radar.

[0011] Optionally, the interferograms are sorted according to the magnitude of the atmospheric phase in each SAR acquisition, starting from the date with the highest noise level. The atmospheric delay before that date is used to correct the interferograms after that date, specifically including: For the temporal edge data of the SAR dataset, the atmospheric delay is estimated and corrected using a spatiotemporal filtering method. The spatiotemporal filtering method is based on the assumption that atmospheric delay is spatially correlated and temporally uncorrelated, and separates the atmospheric delay component in the SAR dataset through spatial low-pass filtering and temporal high-pass filtering. For non-temporal edge data in the SAR dataset, the atmospheric delay is estimated and corrected using the phase stacking method. The phase stacking method takes advantage of the fact that interferograms of the same date contain the same atmospheric delay component, and estimates the atmospheric delay component of the same date by stacking multiple pairs of interferograms.

[0012] This invention provides an atmospheric delay correction system for interferometric synthetic aperture radar, comprising: The data acquisition module is used to acquire Sentinel-1 SAR data captured by interferometric synthetic aperture radar and arrange them according to the capture time sequence to obtain a preprocessed SAR image set. The signal separation module is used to perform time series analysis on the preprocessed SAR image set using StaMPS to obtain multiple interferograms; and to unwrap the multiple interferograms to obtain multiple unwrapped phases; The matrix decomposition module is used to assemble multiple unwrapped phases into an observation matrix. The observation matrix is ​​decomposed by variational Bayesian independent component analysis to distinguish the components representing the layered delay from the independent components by utilizing the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components. The vertical stratification delay correction module is used to identify independent components related to vertical stratification delay based on the correlation coefficient and periodicity coefficient between each independent component and the terrain elevation, and to remove the independent components related to vertical stratification delay from the observation matrix to obtain the residual phase matrix. The turbulence delay correction module is used to perform turbulence delay correction on the residual phase matrix using phase stacking and spatiotemporal filtering methods, so as to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar.

[0013] The atmospheric delay correction method and system for interferometric synthetic aperture radar provided in this invention have the following advantages compared with the prior art: This invention employs variational Bayesian independent component analysis to adaptively decompose the observation matrix composed of unwrapped phases. Instead of presupposing a fixed functional relationship between atmospheric delay and terrain, it utilizes the spatial characteristics (correlated with terrain height) and temporal characteristics (seasonal periodic variations) of vertical stratification delay. By accurately identifying and removing independent components related to vertical stratification delay through the correlation coefficients and periodicity coefficients of each independent component with terrain elevation, it effectively suppresses the interference of atmospheric delay error on deformation signals, preserves the true deformation information to the greatest extent, and significantly improves the accuracy of InSAR deformation monitoring. Attached Figure Description

[0014] Figure 1 This is a flowchart illustrating an atmospheric delay correction method for an interferometric synthetic aperture radar provided in one embodiment; Figure 2 This is an improved CSS schematic diagram of an atmospheric delay correction method for interferometric synthetic aperture radar provided in one embodiment; Figure 3 This is a comparison diagram of vbICA and FastICA decomposition of an atmospheric delay correction method for interferometric synthetic aperture radar provided in one embodiment; Figure 4 The image shows a comparison of phase values ​​before and after interferometric pair correction for an atmospheric delay correction method for an interferometric synthetic aperture radar provided in one embodiment, dated 20210219-20200712. Figure 4 (a) in the image is the original phase map. Figure 4 (b) in the figure is the phase diagram after GACOS correction. Figure 4 (c) in the figure is the phase diagram after linear correction. Figure 4 (d) in the diagram is the phase diagram after ERA5 correction. Figure 4 (e) in the diagram is the phase diagram after FastICA correction. Figure 4 (f) in the figure represents the phase diagram after VBICA correction; Figure 5 The atmospheric delay correction method for interferometric synthetic aperture radar provided in one embodiment is dated 20210219-20220708. A comparison diagram of phase values ​​before and after interferometric pair correction is shown. Figure 5 (a) in the image is the original phase map. Figure 5 (b) in the figure is the phase diagram after GACOS correction. Figure 5 (c) in the figure is the phase diagram after linear correction. Figure 5 (d) in the diagram is the phase diagram after ERA5 correction. Figure 5 (e) in the diagram is the phase diagram after FastICA correction. Figure 5 (f) in the figure represents the phase diagram after VBICA correction; Figure 6 Comparison of residual phase RMS after atmospheric delay correction using different methods in one embodiment of an interferometric synthetic aperture radar; Figure 7 This is a comparison of time series data after atmospheric delay correction using different methods, as shown in one embodiment of an interferometric synthetic aperture radar. Figure 7 (a) in the figure is a comparison of the time series data at location P1 after correction by different methods. Figure 7 (b) in the figure is a comparison of the time series data at location P2 after correction by different methods. Figure 7 (c) in the figure is a comparison of the time series data corrected by different methods at location P3. Figure 7 (d) in the figure is a comparison of the time series data corrected by different methods at location P4. Figure 7 (e) in the figure is a comparison of time series data corrected by different methods at location P5. Figure 7 (f) in the figure is a comparison of time series data corrected by different methods at location P6. Detailed Implementation

[0015] 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.

[0016] In fact, atmospheric errors are mainly influenced by the troposphere and ionosphere. The ionospheric effect is proportional to the signal frequency. Under similar conditions, the ionosphere's influence on the C-band Sentinel-1 satellite is approximately 5.5% of its influence on the L-band satellite. For shortwave SAR data, tropospheric delay is more common in atmospheric errors.

[0017] In existing technologies, the first approach is based on interferometric phase correction. For example, spatiotemporal filtering is a traditional method based on the assumption that the atmosphere is spatially correlated but temporally uncorrelated. However, this assumption is not always valid and easily underestimates the deformed signal. Another approach uses empirical models that utilize the relationship between phase and elevation, which can be described as linear or power-law models. However, it is often difficult to use these models with fixed parameters to mitigate the vertical stratification delay across the entire interferogram. This is because, in most areas with significant elevation variations, the relationship between topography and vertical stratification delay cannot be simply considered linear or power-law.

[0018] The second type is correction based on external data, such as numerical weather models, GNSS, or spectrometer measurements. These external data methods have proven effective in certain specific regions. However, due to the resolution differences between SAR data and weather models, the effectiveness of the correction is often difficult to guarantee.

[0019] Blind Source Separation (BSS), a signal processing and statistical analysis method for solving mixed signal separation problems, has demonstrated its ability to process InSAR time series data and perform atmospheric delay correction. Compared to the traditional independent component analysis algorithm FastICA, Variational Bayesian Independent Component Analysis (vbICA) offers a more flexible method for independent component analysis and exhibits stronger robustness to noise. vbICA allows for adaptive estimation of the probability density distribution of independent components based on variational Bayesian inference. To date, no studies have investigated the application of vbICA separation in terrain-dependent atmospheric delay in interferometric phase analysis.

[0020] Furthermore, after removing vertical stratification delay, the residual tropospheric phase is primarily affected by turbulence. Common Scene Stacking (CSS) is an algorithm designed to estimate turbulent delay without relying on external data. It is based on the premise that interferometric phases with the same date also contain the same atmospheric delay component, and uses a stacking method to estimate the atmospheric delay for the common date. However, due to insufficient constraints at the edges of SAR time series, its estimation of atmospheric delay has high uncertainty. The mechanisms by which atmospheric delay propagation from these edge SAR acquisitions and its impact on deformed time series remain unclear.

[0021] This invention provides an atmospheric delay correction method for interferometric synthetic aperture radar, such as... Figure 1 As shown, the method includes: The Sentinel-1 SAR data captured by interferometric synthetic aperture radar is acquired and arranged according to the capture time sequence to obtain a preprocessed SAR image set.

[0022] StaMPS was used to perform time series analysis on the preprocessed SAR image set to obtain multiple interferograms; and phase unwrapping was performed on the multiple interferograms to obtain multiple unwrapped phases.

[0023] The observation matrix is ​​composed of multiple unwrapped phases. The observation matrix is ​​decomposed by variational Bayesian independent component analysis to distinguish between the components representing the layered delay and the independent components by utilizing the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components.

[0024] Independent components associated with vertical stratification delay are identified based on the correlation coefficients and periodicity coefficients between each independent component and the terrain elevation. The independent components associated with vertical stratification delay are then removed from the observation matrix to obtain the residual phase matrix.

[0025] Turbulent delay correction was performed on the residual phase matrix using phase stacking and spatiotemporal filtering methods to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar.

[0026] A specific embodiment of the present invention is provided: S1: Acquire time-series SAR images of the study area, perform SAR image preprocessing and StaMPS processing (Stanford Method for PS).

[0027] To verify the authenticity and effectiveness of this invention, 39 Sentinel-1 SAR images covering a certain city, collected between January 8, 2020 and December 11, 2022, with a pixel count of 930×5010, were used. After acquiring the time-series SAR images of the monitored area, preprocessing steps such as registration, cropping, and geocoding were performed. Then, StaMPS was used to complete the InSAR time series analysis (StaMPS is a method for performing time-series InSAR processing on preprocessed SAR images; the StaMPS workflow can obtain the unwrapped phase of all interferograms, and the unwrapped phase of all interferogram pairs can be transformed to obtain the time series of deformation in the region, i.e., the deformation situation). The interferogram phase of a pixel is a linear combination of the propagation phase differences from various temporal and spatial sources of variation. (1) in, To untangle the phase, For the deformed phase component, The orbital error phase can be eliminated through linear or quadratic models. This is a residual topographic phase (due to inaccurate DEM). It is random noise. This is the scattering phase (generated due to spatiotemporal decorrelation). This refers to the tropospheric phase.

[0028] For InSAR data, each pixel in the interferogram represents a combination of various noise and deformation sources that vary over time. Independent component analysis (ICA) is a method for separating the mixed signal into its original, statistically independent components. Given that these components are linearly mixed, ICA is well-suited for efficiently decomposing InSAR signals.

[0029] S2: Perform vbICA decomposition based on the unwrapped phase.

[0030] InSAR phase data is generated by t Interferogram and n Composed of individual pixels, it can be organized into an observation matrix. XSubsequently, vbICA decomposition is performed on the interferometric unwrapped phase. The decomposition steps are as follows: vbICA will X Decomposed into a and S A linear combination of these, with added noise: (2) in, For the observation matrix, A It is a mixed matrix. S For the reconstructed source matrix, N It is assumed to be zero-mean Gaussian noise.

[0031] The vbICA model can be described as a model containing a series of parameters, which can be written as... ={ A , Λ , S , q , θ} It is a collection of variables and parameters, where θ ={ Π , μ , β} q It is a discrete variable used to indicate which Gaussian component is contributing to the source signal at a specific location (time point or space point). Λ This is the precision of Gaussian noise, the reciprocal of its variance. μ The center position or average value of each Gaussian component constitutes the source signal distribution. β It is the reciprocal of the variance of each Gaussian component, which describes the shape and dispersion of a single Gaussian component in MoG. Π The weight or probability of each Gaussian component determines the relative contribution of each component to the final source signal distribution. To implement vbICA, it is necessary to calculate the weight or probability of each Gaussian component given the observation matrix. Posterior probability: (3) in, This represents the true posterior probability. As a prior distribution, As the normalization factor, For the variational approximation posterior of each parameter, for The first in One variable, For joint probability, for X The marginal likelihood logarithm, F negative free energy ， for The entropy.

[0032] It is the observation matrix X and W The joint probability is used to quantify the degree of fit between the data and the model. Among them, It can be written as: (4) Since this posterior probability density function is difficult to integrate, variational Bayesian inference is introduced to address this problem. The main concept is to introduce a series of approximate distributions over the weight set. Therefore, This is called the variational distribution, and the Kullback-Leibler (KL) method is used to measure the difference between the true posterior distribution and the variational distribution. (5) (6) in, yes X The marginal likelihood logarithm, F The negative free energy (NFE) is also the objective function. The goal is to maximize the NFE in order to minimize the KL divergence between the true posterior distribution and the variational distribution.

[0033] (7) Next, we can calculate... The optimal solution for the parameters is obtained by taking the partial derivatives with respect to each parameter. By iteratively optimizing the parameters using the Expectation-Maximization (EM) algorithm, the posterior distribution of other unknown sources can be derived.

[0034] because Involves complex integrals, leading to It is difficult to calculate in practice, so variational distribution is introduced. It is an easily computed approximate probability distribution used to approximate the true posterior.

[0035] because Unknown, cannot be minimized directly and KL divergence between Therefore, variational inference is achieved by maximizing the negative free energy. F To indirectly achieve the goal.

[0036] Finally, it can be calculated F The optimal solution for the parameters can be obtained by taking the partial derivative with respect to each parameter. This process is usually iterative until convergence. Finally, the optimal solution is obtained using the mixture matrix. A and source signal S The independent components required for reconstruction.

[0037] The following section will introduce a method for identifying independent components that represent vertical stratification delay.

[0038] S3: Identification of independent components with vertical stratification delay correlation.

[0039] VBICA is used to decompose the phase data. By utilizing the temporal and spatial characteristics of the stratified delay, the components representing the stratified delay are distinguished from other independent components. This method calculates the correlation between independent components and topography, as well as their periodicity in the spectrum, to identify components representing tropospheric delay.

[0040] First, calculate the correlation between each independent component and the terrain. The calculation formula is as follows: (8) in, For the first k The correlation coefficient of each independent component For the first k The first of the independent components i Spatial relationship of pixels, For the first k The mean of the spatial relationships of all pixels of each independent component. H i For the first i The elevation estimated by each InSAR pixel This is the average elevation of all pixels.

[0041] Since this vertical stratification is related to the terrain height, its spatial relationship should theoretically be related to the terrain elevation.

[0042] In addition, the vertical stratification exhibits significant seasonal cyclical variations; therefore, the periodicity of independent components can be estimated by calculating the spectrum. (9) in, For the first k Power spectrum of time series of individual components F For Fourier transform, For the first k The time series of 3 independent components, where ()* represents complex conjugation. L This represents the length of the time series vector.

[0043] Therefore, the spectral features that can be utilized, taking advantage of their highest spatial correlation ( The vertical delay component is determined from the independent components and has strong seasonality. It is then subtracted from the original interferometric synthetic aperture radar LOS measurement to remove the vertical delay component.

[0044] S4: Improved implementation of phase stacking CSS method.

[0045] The CSS method is based on the assumption that there exists a pair of interferograms with the same time interval. For example, consider time... t 1. t 2 and t 3. Two unwrapped interferograms formed by SAR acquisition, with differential phase Δ The formulas associated with each interferogram are as follows: (10) in, It is the date t i Atmospheric delay, Represents the phase of ground deformation, while This represents errors such as those caused by digital elevation models. If we assume that the interferograms have equal time intervals and that deformation occurs at a slow and constant rate, then the differential phase between them is independent of ground deformation. In this case, the same date can be estimated by differentiating the two interferograms. Atmospheric delay, i.e. Δ 12 -Δ 23 .

[0046] If atmospheric phase is uncorrelated with time, combining multiple interferograms from the same date can improve the estimation. The accuracy. With N With the increase in the number of data points, the accuracy of atmospheric phase estimation for each date will improve. (Each data collection date...) The atmospheric delay can be calculated using the following formula: (11) Furthermore, the atmospheric delay for each SAR date needs to be estimated iteratively based on the atmospheric noise figure (ANC). This process first determines the atmospheric delay (APS) for the date with the greatest impact, then removes it to estimate the APS for subsequent dates. Once the estimation is complete, the derived atmospheric phase is subtracted from the interferogram. The ANC is calculated as follows: (12) = (13) in, and This refers to the time sequence number of the SAR image. N The number of SAR image pairs with the same time interval. For the first i Using the background as a reference, the phase difference between consecutive image pairs constructed under the same time interval conditions. For the first j Using the background as a reference, the phase difference between consecutive image pairs constructed under the same time interval conditions. For the first i Scene Pixel m Atmospheric phase delay at that location The average atmospheric phase for all pixels. M This represents the total number of pixels. This represents the RMS value of the atmospheric phase with the highest noise level.

[0047] ANC quantifies the relative level of atmospheric noise in each SAR acquisition and sorts the interferograms according to the magnitude of their atmospheric phase. Starting from the date with the highest noise level, the interferograms of subsequent dates are corrected using the previously estimated atmospheric delay, resulting in the corrected phase and deformation time series of the interferometric synthetic aperture radar.

[0048] However, the accuracy of atmospheric phase estimation often decreases at the edges of SAR acquisition. Therefore, in this invention, CSS is combined with spatiotemporal filtering methods to improve the CSS problem, as follows: Figure 2 As shown, for the temporal edge data of the SAR dataset, a spatiotemporal filtering method is used to estimate atmospheric delay for correction, while for intermediate time data (non-temporal edge data), the CSS method is used for correction.

[0049] Specifically, for the temporal margin data of the SAR dataset, a spatiotemporal filtering method is used to estimate atmospheric delay for correction. This method is based on the assumption that atmospheric delay is spatially correlated and temporally uncorrelated, and separates the atmospheric delay component in the SAR dataset through spatial low-pass filtering and temporal high-pass filtering.

[0050] For non-temporal edge data in the SAR dataset, atmospheric delay is estimated and corrected using the phase stacking method. The phase stacking method leverages the property that interferograms from the same date contain the same atmospheric delay component, estimating the atmospheric delay component for the same date by stacking multiple pairs of interferograms.

[0051] Figures 4 to 6This paper presents a comparison of the interferometric phase decomposition results between the vbICA and FastICA algorithms, and provides time series data of the independent components of tropospheric delay extracted by vbICA. The independent components extracted by vbICA exhibit a stronger and clearer correlation with the terrain. In contrast, the tropospheric delay signal identified by FastICA is scattered across multiple independent components, leading to a more complex physical interpretation. For ease of visualization, the time series data has been normalized and presented along with its spectrum. The results show a clear annual periodicity, which is highly consistent with the typical characteristics of vertical stratification delay, further validating the advantage of vbICA over FastICA in extracting and isolating vertical stratification delay signals.

[0052] Figure 4 Table 1 further validates the effectiveness of vbICA in vertical layer delay correction. Figure 5 and Figure 6 The residual phase distributions of some interferograms are shown in their original state and after correction using various methods. Theoretically, in unvarnished interferometric phases, the closer the corrected phase is to zero, the better the correction effect. It can be observed that the vbICA method significantly removes more vertical layered delay signals, while other methods all exhibit varying degrees of residual phase.

[0053] To quantify the overall corrective effect of each method, Figure 7 Table 1 presents a quantitative comparison of the root mean square (RMS) of the residual phase in all interferograms and the percentage improvement in the interferograms. A lower RMS indicates a better correction effect. The results show that the present invention achieves the lowest RMS in most interferograms, with an average correction rate of 47.86%, and 94.74% of the interferograms show improved phase after correction, which is significantly better than other methods.

[0054] Table 1. Quantitative comparison of residual phase and percentage improvement in interferogram after correction by different methods.

[0055] Furthermore, to further evaluate the performance of the proposed method in deformation monitoring, typical monitoring points were selected to compare the time-series deformation results under different methods. For example... ​ As shown, the original time series exhibits significant seasonal cyclical variations due to the influence of vertical stratification delay. The correction effects of different methods vary considerably. This invention effectively removes seasonal variations and residual random turbulence errors to the greatest extent possible. Furthermore, compared to the original CSS method, the improved CSS method results in a more stable time series after correction, effectively mitigating the problem of inaccurate turbulence delay estimation at both ends of the time series.

[0056] Based on the same inventive concept, embodiments of the present invention provide an atmospheric delay correction system for interferometric synthetic aperture radar, the system comprising: The data acquisition module is used to acquire Sentinel-1 SAR data captured by interferometric synthetic aperture radar and arrange them according to the capture time sequence to obtain a preprocessed SAR image set.

[0057] The signal separation module is used to perform time series analysis on the preprocessed SAR image set using StaMPS to obtain multiple interferograms; and to unwrap the multiple interferograms to obtain multiple unwrapped phases.

[0058] The matrix decomposition module is used to assemble multiple unwrapped phases into an observation matrix. The observation matrix is ​​decomposed by variational Bayesian independent component analysis to distinguish between components representing the layered delay and independent components using the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components.

[0059] The vertical stratification delay correction module is used to identify independent components related to vertical stratification delay based on the correlation coefficient and periodicity coefficient between each independent component and the terrain elevation, and remove the independent components related to vertical stratification delay from the observation matrix to obtain the residual phase matrix.

[0060] The turbulence delay correction module is used to perform turbulence delay correction on the residual phase matrix using phase stacking and spatiotemporal filtering methods, so as to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar.

[0061] The embodiments described above are merely examples of several implementations of the present invention, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. An atmospheric delay correction method for interferometric synthetic aperture radar, characterized in that, include: Acquire Sentinel-1 SAR data captured by interferometric synthetic aperture radar and arrange them according to the capture time sequence to obtain a preprocessed SAR image set; StaMPS was used to perform time series analysis on the preprocessed SAR image set to obtain multiple interferograms; and phase unwrapping was performed on the multiple interferograms to obtain multiple unwrapped phases; The observation matrix is ​​composed of multiple unwrapped phases. The observation matrix is ​​decomposed by variational Bayesian independent component analysis to distinguish the components representing the layered delay from the independent components by utilizing the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components. The independent components related to vertical stratification delay are identified based on the correlation coefficient and periodicity coefficient between each independent component and the terrain elevation. The independent components related to vertical stratification delay are then removed from the observation matrix to obtain the residual phase matrix. Turbulent delay correction was performed on the residual phase matrix using phase stacking and spatiotemporal filtering methods to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar.

2. The atmospheric delay correction method for interferometric synthetic aperture radar as described in claim 1, characterized in that, The process of unwrapping multiple interferograms to obtain multiple unwrapped phases specifically includes: Each interferogram is considered as a linear combination of the propagating phase differences between various temporal and spatial variation sources, based on the following formula: ； in, To untangle the phase, For the deformed phase component, For orbital error phase, For residual topographic phase, It is random noise. For the scattering phase, This refers to the tropospheric phase.

3. The atmospheric delay correction method for interferometric synthetic aperture radar as described in claim 1, characterized in that, The observation matrix is ​​composed of multiple unwrapped phases. Variational Bayesian independent component analysis is then used to decompose the observation matrix to distinguish between components representing the hierarchical delay and independent components using the temporal and spatial characteristics of the hierarchical delay, resulting in multiple independent components, specifically including: The observation matrix is ​​decomposed using variational Bayesian independent component analysis based on the following formula: ； in, For the observation matrix, A It is a mixed matrix. S For the reconstructed source matrix, N Assuming it is zero-mean Gaussian noise; The mixing matrix and the reconstructed source matrix are described as the vbICA model based on the following formula: ={ A , Λ , S , q , θ}； θ ={ Π ， μ ， β}； in, A collection of variables and parameters. Λ For Gaussian noise accuracy, q For discrete variables, μ The center location or average value of each Gaussian component. β The reciprocal of the variance of each Gaussian component. Π The weight or probability of each Gaussian component; Given the observation matrix, calculate Posterior probability: ； ； ； ； ； in, This represents the true posterior probability. As a prior distribution, As the normalization factor, For the variational approximation posterior of each parameter, for The first in One variable, For joint probability, for X The marginal likelihood logarithm, F negative free energy ， for Entropy; Multiple independent components are constructed based on the mixing matrix and the reconstructed source matrix.

4. The atmospheric delay correction method for interferometric synthetic aperture radar as described in claim 1, characterized in that, The correlation coefficients between each independent component and topographic elevation are determined based on the following formula: ； in, For the first k The correlation coefficient of each independent component For the first k The first of the independent components i Spatial relationship of pixels, For the first k The mean of the spatial relationships of all pixels of each independent component. H i For the first i The elevation estimated by each InSAR pixel This is the average elevation of all pixels; The periodicity coefficients between each independent component and the terrain elevation are determined based on the following formula: ； in, For the first k Power spectrum of time series of individual components F For Fourier transform, For the first k The time series of 3 independent components, where ()* represents complex conjugation. L This represents the length of the time series vector.

5. The atmospheric delay correction method for interferometric synthetic aperture radar as described in claim 1, characterized in that, The process of correcting the residual phase matrix for turbulence delay using phase stacking and spatiotemporal filtering methods to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar specifically includes: The atmospheric delay for each acquisition date in the residual phase matrix is ​​estimated using the phase stacking method based on the following equation: ； The atmospheric delay for each SAR date is estimated iteratively based on the atmospheric noise figure using the following formula: ； = ； in, and This refers to the time sequence number of the SAR image. N The number of SAR image pairs with the same time interval. For the first i Using the background as a reference, the phase difference between consecutive image pairs constructed under the same time interval conditions. For the first j Using the background as a reference, the phase difference between consecutive image pairs constructed under the same time interval conditions. For the first i Scene Pixel m Atmospheric phase delay at that location The average atmospheric phase for all pixels. M This represents the total number of pixels. The RMS value of the atmospheric phase with the highest noise level; The interferograms are sorted according to the magnitude of the atmospheric phase in each SAR acquisition. Starting from the date with the highest noise level, the atmospheric delay before that date is used to correct the interferograms after that date, resulting in the phase and deformation time series corrected by the interferometric synthetic aperture radar.

6. The atmospheric delay correction method for interferometric synthetic aperture radar as described in claim 5, characterized in that, The process involves sorting the interferograms according to the magnitude of the atmospheric phase in each SAR acquisition, starting from the date with the highest noise level, and using the atmospheric delay before that date to correct the interferograms after that date. Specifically, this includes: For the temporal edge data of the SAR dataset, a spatiotemporal filtering method is used to estimate atmospheric delay for correction. The spatiotemporal filtering method is based on the assumption that atmospheric delay is spatially correlated and temporally uncorrelated, and separates the atmospheric delay component in the SAR dataset through spatial low-pass filtering and temporal high-pass filtering. For non-temporal edge data in the SAR dataset, atmospheric delay is estimated and corrected using a phase stacking method. This phase stacking method utilizes the characteristic that interferograms from the same date contain the same atmospheric delay component, and estimates the atmospheric delay component from the same date by stacking multiple pairs of interferograms.

7. An atmospheric delay correction system for interferometric synthetic aperture radar, characterized in that, include: The data acquisition module is used to acquire Sentinel-1 SAR data captured by interferometric synthetic aperture radar and arrange them according to the capture time sequence to obtain a preprocessed SAR image set. The signal separation module is used to perform time series analysis on the preprocessed SAR image set using StaMPS to obtain multiple interferograms; and to unwrap the multiple interferograms to obtain multiple unwrapped phases; The matrix decomposition module is used to assemble multiple unwrapped phases into an observation matrix. The observation matrix is ​​decomposed by variational Bayesian independent component analysis to distinguish the components representing the layered delay from the independent components by utilizing the temporal and spatial characteristics of the layered delay, thus obtaining multiple independent components. The vertical stratification delay correction module is used to identify independent components related to vertical stratification delay based on the correlation coefficient and periodicity coefficient between each independent component and the terrain elevation, and to remove the independent components related to vertical stratification delay from the observation matrix to obtain the residual phase matrix. The turbulence delay correction module is used to perform turbulence delay correction on the residual phase matrix using phase stacking and spatiotemporal filtering methods, so as to obtain the corrected interferometric phase and deformation time series of the interferometric synthetic aperture radar.

Images