An InSAR atmospheric delay correction method and system based on adaptive window

Abstract

The present invention discloses an InSAR atmospheric delay correction method and system based on an adaptive window. The method includes the following steps: step S1, data acquisition and preprocessing; step S2, InSAR data processing; step S3, adaptive dynamic window segmentation; step S4, multi-temporal stratospheric delay, deformation and terrain residual joint modeling; step S5, model parameter M estimation; step S6, unwrapped phase correction and cumulative deformation sequence solution. A system for an InSAR atmospheric delay correction method based on an adaptive window. The present invention introduces an atmospheric moisture content product and performs adaptive SAR image window segmentation based on the difference in water vapor content in the spatial atmosphere. It also considers the relationship between vegetation cover density and interferogram coherence, and uses the average coherence of the InSAR interferogram as the initial weight matrix, thereby reducing the impact of seasonal vegetation density changes on deformation results. The M estimation method is used to iteratively solve the InSAR deformation time series, thereby improving the reliability of the deformation estimation results.

Description

Technical Field

[0001] The present invention relates to an InSAR atmospheric delay correction method and system, and in particular to an InSAR atmospheric delay correction method and system based on an adaptive window. Background Art

[0002] InSAR (Interferometric Synthetic Aperture Radar) is a new space-based geodetic and remote sensing tool. Combining synthetic aperture radar (SAR) and interferometry, it offers advantages such as all-day, all-weather, and continuous, large-scale observations. It has been widely used to monitor natural disasters such as earthquakes, landslides, and glacial movement, as well as human activities. However, the demands of economic development, such as investigating ground subsidence hazards along high-speed rail lines and monitoring landslides and dangerous rock masses, have placed higher demands on InSAR monitoring accuracy. However, atmospheric delay, a common tropospheric delay, can result in deformation errors of tens of centimeters, hindering further improvements in InSAR monitoring accuracy. Therefore, research on how to mitigate or correct the effects of atmospheric delay in InSAR is crucial to ensuring the accuracy of monitoring results. This not only helps fully tap the potential of InSAR technology in fields such as geological hazard investigation but also has important implications for ensuring the safety and reliability of economic development projects.

[0003] Currently, most InSAR atmospheric delay phase correction methods employ corrections based on phase-elevation models. For example, these methods employ a fixed-size sliding window to construct a model of the relationship between tropospheric delay and terrain elevation within a local area, using elevation values ​​to simulate tropospheric delay. These methods often rely on empirical judgment to segment SAR (Synthetic Aperture Radar) images, which compromises the accuracy of the atmospheric delay model. Furthermore, improper SAR image segmentation can lead to computational redundancy.

[0004] In summary, the current nSAR atmospheric delay method has the following main shortcomings: ① Traditional local window segmentation usually adopts a fixed size and does not consider the heterogeneity of atmospheric spatial distribution, which affects the accuracy of the atmospheric delay model and causes computational redundancy; ② SAR satellites mainly use shortwave bands, which have relatively weak penetration ability into vegetation. Seasonal changes in vegetation density can affect the coherence of the interferogram in InSAR data processing, which can lead to more errors or local data missing in the interferogram, affecting the accuracy of the deformation results. Summary of the Invention

[0005] In order to address the shortcomings of the above-mentioned technologies, and to address the problem that phase delays caused by differences in the atmospheric state on temporal and spatial scales easily mask the true deformation of the surface, the present invention provides an InSAR atmospheric delay correction method and system based on an adaptive window, which can perform atmospheric delay correction on InSAR data, thereby ensuring the accuracy of the time series deformation results.

[0006] In order to solve the above technical problems, the technical solution adopted by the present invention is: an InSAR atmospheric delay correction method based on an adaptive window, the method comprising the following steps:

[0007] Step S1, data acquisition and preprocessing;

[0008] Step S2: InSAR data processing:

[0009] Step S3, adaptive dynamic window segmentation: using the variogram to calculate the semivariogram value list, the number of SAR image segmentation windows is determined by the size of the model range obtained by fitting, and the SAR image window segmentation is completed using the K-means clustering method;

[0010] Step S4, joint modeling of multi-temporal stratum delay, deformation and terrain residual;

[0011] Step S5, estimating model parameters M;

[0012] Step S6, unwrapping phase correction and cumulative deformation sequence solution: the average coherence of the interferogram is introduced as the initial weight matrix, and the M estimation method is used again to solve the deformation time series to achieve atmospheric delay correction of the deformation results.

[0013] Preferably, in step S1, the following data need to be acquired: multi-temporal InSAR data, ERA5 data corresponding to the time when the InSAR data was acquired, and DEM data; the acquired data are preprocessed:

[0014] Multi-temporal InSAR data preprocessing: including SLC extraction, DEM geocoding, image registration, study area cropping, and multi-viewing;

[0015] ERA5 data preprocessing: This includes data format conversion, spatial region cropping, vertical integration of temperature and relative humidity to obtain the atmospheric water vapor content of the entire troposphere, and then calculating the average atmospheric water vapor content within the spatial region. Spatial interpolation is also performed for cases where the spatial resolution is insufficient.

[0016] Preferably, in step S2, the pre-processed SLC image is first subjected to interferometric baseline calculation according to the time baseline threshold and the spatial baseline threshold. Then, based on the connection of the baselines, the SAR image is subjected to differential interferometry, filtering, phase unwrapping, and elevation correction. The phase formula after unwrapping is as follows:

[0017] Φ=Φ defo +Φ topo +Φ strat +Δ (1)

[0018] Where Φ represents the unwrapped phase of the entire interferogram, Φ defo represents the deformation phase, Φ topo represents the terrain residual phase, Φ strat represents the tropospheric delay phase, and Δ represents the error vector.

[0019] Preferably, step S3 specifically includes the following steps:

[0020] Step 1: Calculation of variance function

[0021] According to the spatial distribution of atmospheric moisture, the semivariogram function value list of atmospheric water vapor content and spatial position data is calculated using the variogram, and the spherical model or exponential model is selected for fitting to obtain the model range;

[0022] Step 2: Window division and size adjustment

[0023] The number of SAR image segmentation windows is determined according to the range size. The longitude, latitude and atmospheric water vapor content are classified using the K-means clustering method to obtain multiple clusters of different types. The longitude and latitude coordinates at the boundary of each cluster are obtained to complete the SAR image window segmentation.

[0024] Preferably, step one specifically includes the following steps:

[0025] Based on the average atmospheric water vapor content and spatial location data obtained from ERA5 data, the semivariogram function is used as the variation function to describe the degree of change in water vapor content between different locations. The specific calculation formula is:

[0026]

[0027] Among them, χ(d) is the semivariogram function value, z(x i ) and z(x j ) are respectively at the spatial position x i 、x j The atmospheric water vapor content at , d is the spatial distance, and N(d) is the total number of pixel pairs with a spatial distance of d;

[0028] Use the Euclidean distance formula to calculate the distance between each pixel and its coordinates (xi ,y i ) and (x j ,y j ) Calculate the spatial distance d:

[0029]

[0030] According to the calculated spatial distance d, all pixel pairs are classified and the pixel pairs within the same distance range are grouped into one category. For each category, i.e., pixel pairs with the same d value, the semivariance function value is calculated according to formula (2). The semivariance function value corresponding to each distance d is recorded to form a semivariance function value list.

[0031] Preferably, step 2 specifically includes the following steps:

[0032] Use the semivariogram value list calculated in step 1 to determine the window size, select an appropriate semivariogram model for fitting, and obtain the model range.

[0033] According to the range value and the range of the SAR image processed in step S2, that is, the size of the SLC image obtained by data preprocessing, the number of windows in the area is preliminarily determined. Then, the longitude, latitude, and water vapor content values ​​of each pixel in the SAR image are used as the sample set, and the sample set is divided into multiple clusters (U1, U2, ..., U u ), each cluster must satisfy the minimum value of the square error E:

[0034]

[0035] Where α represents the i-th cluster U i The sample value vector corresponding to each pixel in ; represents the i-th cluster U i The mean vector of , often called the centroid, is expressed as U i | is the number of samples in the i-th cluster;

[0036] After the sample set is divided, the centroid and longitude and latitude values ​​at the boundary of each cluster are extracted, and the coordinate vector of the inflection point at the boundary of each window is determined according to the following formula:

[0037]

[0038] in, The new boundary inflection point coordinate vector is calculated to ensure that there is a 20% overlap between the i-th cluster and the adjacent cluster; represents the coordinate vector of the inflection point at the boundary of the i-th cluster that satisfies the conditions of formula (4); Represents the coordinates of the centroid position of the i-th cluster; In summary, the coordinate vector of the boundary inflection point of each cluster The SAR image segmentation can be completed.

[0039] Preferably, step S4 specifically includes the following steps:

[0040] Step 1: Build a multi-time relative spheric delay model window by window

[0041] Based on the window segmented in step S3, a linear model of tropospheric atmospheric delay phase and elevation is constructed pixel by pixel in any single window for a single interferogram:

[0042]

[0043] in, H represents the tropospheric delay phase difference between the nth and n+1th SAR images, hgt is the elevation parameter matrix corresponding to the pixels in the window, Representatives and The corresponding linear model parameter matrix of phase and elevation; s is the total number of pixels in a single window, represents the tropospheric delay phase difference of the qth pixel (1≤q≤s); h q Indicates the elevation value corresponding to the qth pixel (1≤q≤s); k n ,ξ n represent the slope and intercept of the linear model of phase and elevation, respectively, and T represents the matrix or vector transpose;

[0044]

[0045] Assume that InSAR data processing acquires N SAR images covering the same study area and generates M differential interferograms, where M satisfies Assuming that the tropospheric delay of the reference image is zero, the unwrapped phase of M×s pixels in a single window of M interferograms is obtained according to formula (1), which is specifically expressed as follows:

[0046]

[0047] in, It represents the unwrapped phase of M×s pixels in the same window of M interference patterns, A represents the matrix related to the organization of the interference pattern, Indicates the calculation of the Kronecker product of two matrices, H hgt is the elevation parameter matrix corresponding to the s pixels in the window, X strat The linear model parameter matrix representing the tropospheric atmospheric delay phase and elevation of any adjacent SAR images, They represent the deformation phase and topographic residual phase in the M interferograms, which will be modeled in the second step. Δ represents the error vector, which is estimated by residual analysis or using an external error model.

[0048] is the unwrapped phase of the qth pixel in the mth interferogram (1≤m≤M), h q Indicates the elevation value corresponding to the qth pixel (1≤q≤s), k n ,ξ n represent the slope and intercept of the linear model of the tropospheric delay phase and elevation between the nth and n+1th SAR images (1≤n<N);

[0049]

[0050] Step 2: Construct surface deformation and terrain residual models window by window

[0051] A cubic polynomial is used to describe the temporal variation of surface deformation. At the same time, the terrain residual is modeled based on the geometric model of the InSAR terrain phase. Therefore, the specific formulas for the deformation phase and terrain residual phase of each pixel in a single window are as follows:

[0052]

[0053] Among them, t i , t j represent the imaging time of the i-th and j-th SAR images respectively; Respectively represent t i and t j The surface deformation phase and terrain residual phase of the qth pixel in time;

[0054] Δh q are the parameters to be estimated, representing the average velocity, acceleration, acceleration change rate and terrain residual phase of the qth pixel respectively; λ is the radar wavelength; r is the slant range between the SAR sensor and the ground, and θ is the incident angle of the SAR satellite sensor. Both r and θ are obtained from the SLC parameter file; Indicates t i , t j vertical baseline between time sensors;

[0055] Extending formula (8) to M interferograms, the surface deformation phase and terrain residual phase of M×s pixels in a single window are expressed as follows:

[0056]

[0057] Among them, I s×sRepresents the unit matrix of size s×s; the matrix G represents the time matrix corresponding to the matrix A, which is expressed as M represents the total number of interferograms, N represents the total number of SAR images, and t n+1 -t n is the time interval between adjacent SAR images (1≤n<N);

[0058] D defo 、D topo They represent the model coefficient matrices related to the surface deformation phase and the terrain residual phase, respectively. The specific form is expressed as:

[0059]

[0060] Among them, t n+1 -t n is the time interval between adjacent SAR images (1≤n<N), Δh q represents the terrain residual phase of the qth pixel (1≤q≤s), X defo 、X topo Represent the parameter vectors to be estimated respectively;

[0061]

[0062] Step 3: Joint modeling of multi-temporal relativity spheric delay, surface deformation, and terrain residuals

[0063] Substituting formula (9) into formula (7), we can construct the observation equation of M interference patterns with M×s number of pixels in a single window:

[0064]

[0065] in, is the actual observation value, representing the unwrapped phase of M×s pixels, B represents the coefficient matrix of the observation equation, and X represents the parameter vector to be estimated, which is expressed as X=[X start ,X defo ,X topo ] T ;

[0066]

[0067] Preferably, in step S5, the M estimation is used to iteratively solve the parameter vector X to be estimated and the InSAR deformation time series in formula (10), including the following two steps:

[0068] Step 1: Calculate the initial weight matrix

[0069] Calculate the overall mean of the complex signal sum in the local area centered on pixel g. The specific expression is:

[0070]

[0071] Among them, γ m,g represents the coherence of the g-th pixel in the m-th interference pattern (1≤m≤M,1≤g≤T), where T represents the total number of interference image elements; c i (μ) and c j (μ) represents the complex signal of each pixel in the local area centered on pixel g in the i-th and j-th SAR images, respectively. Both are obtained from the SLC file. i and j are the two SAR images used to generate the interferogram m (1≤i<j≤N); R is the number of pixels in the local area centered on pixel g; * is a conjugate complex number;

[0072] After calculating the coherence of each pixel in the interference pattern pixel by pixel using formula (11), the weight of the corresponding interference pattern is determined by the following formula:

[0073]

[0074] Among them, p m Indicates the weight of the mth interference pattern; represents the average coherence of the mth interference pattern, expressed as γ m,g is the coherence corresponding to each pixel calculated by formula (11) (1≤g≤T); C represents the average coherence of the unit weight interference pattern;

[0075] Therefore, the initial weight matrix P is expressed as:

[0076]

[0077] Step 2: M estimation algorithm process

[0078] The solution of the parameter vector X under M estimation is obtained through formula (10). First, the ρ function is introduced as the loss function to optimize the parameters:

[0079]

[0080] Among them, ρ(x) is a real-valued function and is not unique, expressed as Usually ω=1.345; represents the row vector, representing the ii-th row of the coefficient matrix B of the observation equation (1≤ii≤M×s); is the M-estimate solution of the parameter vector X; v ii is the observed value φ ii Correction number; p ii is the weight of the ii-th observation, calculated from formula (13); ii Represents the ii-th observation value, which is a column vector The unwrapped phase value of the value at the ii-th position in ;

[0081] The column vector is represented as

[0082] If the order After derivation of formula (14), the solution equation of M estimation is:

[0083]

[0084] Therefore, in the case of multiple parameters, the parameters Solve it iteratively using the following formula:

[0085]

[0086] Where W is the equivalent weight matrix, which is constructed using the ψ(x) function and the ρ(x) function, specifically expressed as W = diag(W1, W2, ..., W M×s ),in or

[0087] Preferably, in step S6, after completing the parameter estimation of the joint model of the unwrapped phase and tropospheric delay, deformation and terrain residual in a single window, the following steps will be continued:

[0088] S61, window-by-window parameter estimation: executing steps S4 and S5 to complete the establishment and parameter estimation of the tropospheric delay model for all windows determined in step S3 one by one;

[0089] S62, parameter space interpolation and tropospheric delay reconstruction: First, the obtained tropospheric delay parameter values ​​are recorded at the positions of the geometric center pixels of their respective windows. Then, the parameter values ​​at the geometric center positions of each window in each image sheet are used to interpolate the entire image sheet to obtain the parameter values ​​corresponding to all pixel positions in each image sheet. Finally, the tropospheric delay phase is reconstructed according to the elevation value corresponding to each pixel using formula (6).

[0090] S63, converting the tropospheric delay phase: multiplying the reconstructed tropospheric delay phase by a conversion matrix A on the left to convert the tropospheric delay phase into a form corresponding to the interferogram connection mode;

[0091] S64, unwrapped phase correction: removing the converted tropospheric delay phase and terrain residual phase from the unwrapped phase to complete the tropospheric atmospheric delay correction of the unwrapped phase;

[0092] S65. Deformation time series solution: Use the average coherence of the interferogram as the initial weight matrix and apply the M estimation method of step S5 again to the unwrapped phase Φ′ after atmospheric delay correction: Replace the coefficient matrix B in formula (16) with the matrix A, and iteratively solve the cumulative deformation phase in any adjacent SAR imaging period

[0093] A system for InSAR atmospheric delay correction method based on an adaptive window, the system comprising:

[0094] Data preprocessing module: uses ERA5 data to calculate the atmospheric water vapor content and its spatial location;

[0095] Window adaptive dynamic segmentation module: SAR image window segmentation is completed based on the value list obtained by the data preprocessing module according to the variance function calculation and combined with the K-means clustering method;

[0096] InSAR data processing module: The interferometric baseline is calculated for the SLC image obtained by the data preprocessing module according to the set time baseline threshold and spatial baseline. Based on the connection of the baseline, the SAR image is subjected to differential interferometry, filtering, phase unwrapping, and elevation correction to obtain the interferogram unwrapping phase.

[0097] Multi-temporal relative spheric delay, deformation, and terrain residual joint modeling module: This module combines the segmentation window determined by the window adaptive dynamic segmentation module and the interferogram unwrapped phase obtained by the InSAR data processing module, and uses the elevation value, incident angle, slant range, and vertical baseline corresponding to each pixel in the window to construct a multi-temporal relative spheric delay, surface deformation, and terrain residual joint model window by window;

[0098] M estimation algorithm module: This module provides two initial weight matrices. The first is the unit weight matrix. The second is to calculate the coherence between the pixels corresponding to the two SAR images of each interferogram based on the baseline connection information in the InSAR data processing module and the complex signal of each pixel in the SAR image. Finally, the average coherence of each interferogram is calculated as the initial weight matrix. The estimated parameters are solved through continuous iteration to achieve adaptive adjustment of the residual.

[0099] Unwrapped phase tropospheric delay correction module: This module first uses the unit weight matrix of the M estimation algorithm module as the initial weight matrix to iteratively solve the parameter model parameters, then uses the model parameters to reconstruct the tropospheric delay phase and terrain residual phase, and finally removes this part of the phase from the unwrapped phase to complete the unwrapped phase tropospheric delay correction;

[0100] Cumulative deformation time series solution module: The unwrapped phase of the interferogram corrected by the unwrapped phase tropospheric delay correction module is used as the observation value, and the average coherence of the interferogram provided by the M estimation algorithm module is used as the initial weight matrix to iteratively solve the cumulative deformation time series.

[0101] Compared with the prior art, the present invention has the following beneficial effects:

[0102] First, based on the spatial distribution of atmospheric moisture content, the K-Means clustering method is used to adaptively segment SAR images into windows based on longitude, latitude, and spatial atmospheric water vapor content. By modeling each local window, the layered atmospheric delay phase during InSAR processing can be accurately obtained, overcoming the inaccurate modeling problem caused by atmospheric spatial heterogeneity in traditional fixed-size window methods.

[0103] Second, considering the relationship between vegetation density changes and interferogram coherence, the average coherence of the InSAR interferogram is used as the initial weight matrix, and combined with M estimation to solve the InSAR deformation time series. This can compensate for the unreliable deformation results caused by seasonal vegetation density changes, and ensure more accurate estimation results through adaptive adjustment of residuals.

[0104] The present invention targets the spatial distribution differences of atmospheric water vapor content and introduces a variogram to quantitatively describe the spatial heterogeneity of the atmosphere. Based on this, the K-Means clustering method is used to perform adaptive window segmentation based on longitude, latitude, and spatial atmospheric water vapor content. By modeling each local window, the layered atmospheric delay phase can be accurately obtained in InSAR data processing while avoiding computational redundancy, thereby more accurately modeling and correcting the atmospheric layered delay.

[0105] The present invention comprehensively considers the relationship between vegetation cover density and interferogram coherence, uses the average coherence of the InSAR interferogram as the initial weight matrix, and combines the M estimation method to iteratively solve the InSAR deformation time series. It can effectively reduce the impact of vegetation density changes on the deformation results, thereby improving the accuracy of surface deformation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0106] Figure 1 It is the overall flow chart of the present invention.

[0107] Figure 2 It is a system composition diagram of the present invention. DETAILED DESCRIPTION

[0108] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0109] Aiming at the problem that deformation results in InSAR data processing are easily affected by atmospheric delay, the present invention proposes a method and system for InSAR atmospheric delay correction based on adaptive window. First, the spatial heterogeneity of atmospheric moisture content is quantitatively described. According to the spatial distribution of atmospheric moisture, the semivariogram function value list of atmospheric water vapor content and spatial position data is calculated using the variation function, and a spherical model or an exponential model is selected for fitting to obtain the model variation range.

[0110] Then, the number of SAR image segmentation windows is determined according to the range size. The K-means clustering method is used to classify the longitude, latitude and atmospheric water vapor content to obtain multiple clusters of different types. The longitude and latitude coordinates at the boundary of each cluster are obtained to complete the SAR image window segmentation. Then, a relationship model between the InSAR unwrapped phase and the deformation signal, tropospheric delay and terrain residual is established in each local window. The M estimation method is used to solve the parameters in the model with unit weight as the initial weight matrix, and the tropospheric delay within the imaging time of adjacent SAR images is estimated to complete the atmospheric delay correction of the interferogram unwrapped phase.

[0111] Finally, for the corrected unwrapped phase, the average coherence of the interferogram is introduced as the initial weight factor, and the M-estimation method is used again to solve the deformation time series to achieve atmospheric delay correction of the deformation results.

[0112] The overall process of the present invention is as follows Figure 1 As shown in the figure, the main steps are as follows: 1) data acquisition and preprocessing; 2) InSAR data processing; 3) adaptive dynamic window segmentation; 4) joint modeling of multi-temporal relative stream layer delay, deformation and terrain residual; 5) estimation of model parameter M; 6) unwrapped phase correction and cumulative deformation sequence solution.

[0113] (1) Data acquisition and preprocessing

[0114] In order to perform InSAR layered delay correction and deformation sequence solution, the following data are required:

[0115] 1) Multi-temporal InSAR data (such as Sentinel-1A SLC products), used for InSAR data processing experiments;

[0116] 2) ERA5 reanalysis data: It is necessary to select ERA5 data corresponding to the acquisition time of InSAR data. This product will be used to calculate the spatial distribution heterogeneity of the atmosphere to determine the adaptive window;

[0117] 3) DEM data (such as 30m resolution STRM data), which is used to remove the terrain phase of the interference pattern during InSAR data processing.

[0118] The preprocessing of the above data includes the following steps:

[0119] 1) Sentinel-1A SLC product preprocessing: including SLC extraction, DEM geocoding, image registration, study area cropping, and multi-viewing;

[0120] 2) ERA5 reanalysis data preprocessing: This includes data format conversion, spatial region cropping, vertical integration of temperature and relative humidity to obtain the atmospheric water vapor content for the entire troposphere, and calculation of the average atmospheric water vapor content within the spatial region. Spatial interpolation can be performed to address insufficient spatial resolution. Quality control is performed on the preprocessed ERA5 data to remove any abnormal or unreasonable data points.

[0121] These preprocessing steps provide the necessary input data for subsequent InSAR layer delay correction and ensure the quality and applicability of the data.

[0122] (2) InSAR data processing

[0123] After preprocessing, the Sentinel-1A SLC product first calculates an interferometric baseline using a temporal baseline threshold of 36 days and a spatial baseline threshold of 300 meters. These temporal and spatial baseline thresholds can be adjusted based on the SAR data. Then, based on the baseline connectivity, the SAR imagery undergoes differential interferometry, filtering, phase unwrapping, and elevation correction. The unwrapped phase formula is as follows:

[0124] Φ=Φ defo +Φ topo +Φ strat +Δ (1)

[0125] Among them, Φ defo represents the deformation phase, Φ topo represents the terrain residual phase, Φ strat represents the tropospheric delay phase, and Δ represents the error vector. Here, Φ represents the unwrapped phase of the entire interferogram, which will be used as the observation value to solve the model parameters in the subsequent steps.

[0126] (3) Adaptive dynamic window segmentation

[0127] Step 1: Calculation of the variance function

[0128] First, based on the average atmospheric water vapor content and spatial location data obtained from ERA5 data, the semivariogram function is used as the variation function to describe the degree of change in water vapor content between different locations. The specific calculation formula can be expressed as:

[0129]

[0130] Here χ(d) is the semivariogram function value, z(x i ) and z(x j ) are respectively at the spatial position x i 、x j The atmospheric water vapor content at , d is the spatial distance, and N(d) is the total number of pixel pairs with a spatial distance of d.

[0131] For each pixel pair coordinate (x i ,y i ) and (x j ,y j ), to calculate the spatial distance d, the Euclidean distance formula can be used:

[0132]

[0133] According to the calculated spatial distance d, all pixel pairs are classified and the pixel pairs within the same distance range are grouped into one category.

[0134] For each class (i.e., pixel pairs with the same d value), calculate its semivariance function value according to formula (2), record the semivariance function value corresponding to each distance d, and form a semivariance function value list.

[0135] Step 2: Window division and size adjustment

[0136] The window size is determined using the list of semivariogram values ​​calculated in the first step. An appropriate semivariogram model, such as an exponential or spherical model, is selected for fitting. Model parameters, including the sill, nugget, and range, can be obtained using fitting methods such as nonlinear least squares. The fitted model range represents the distance from the starting point to where the semivariogram reaches a stable value, reflecting the range of spatial correlation—that is, the range within which the data are correlated.

[0137] According to the range value, combined with the range of the SAR image processed in step (2), that is, the size of the SLC image obtained by data preprocessing (the SAR image range and the SLC image range are consistent here. In the InSAR data processing process, the SLC image is used as input data. Therefore, the range of the SAR image involved in the subsequent InSAR processing process is the same as the size of the SLC image), the number of windows in the area can be preliminarily determined. Then, the longitude, latitude, and water vapor content values ​​of each pixel in the SAR image are used as the sample set. According to the distance between samples, the sample set can be divided into multiple clusters (U1, U2, ..., U u ), each cluster must satisfy the minimum value of the square error E:

[0138]

[0139] Where α represents the i-th cluster Ui The sample value vector corresponding to each pixel in ; represents the i-th cluster U i The mean vector of , often called the centroid, can be expressed as U i | is the number of samples in the i-th cluster.

[0140] After completing the division of the sample set, extract the centroid and longitude and latitude values ​​at the boundary of each cluster, and determine the coordinate vector of the inflection point at the boundary of each window according to the following:

[0141]

[0142] here The new boundary inflection point coordinate vector is calculated to ensure that there is a 20% overlap between the i-th cluster and the adjacent cluster; represents the coordinate vector of the inflection point at the boundary of the i-th cluster that satisfies the conditions of formula (4); Represents the coordinates of the centroid position of the i-th cluster; In summary, the coordinate vector of the boundary inflection point of each cluster That is, SAR image segmentation can be completed.

[0143] (4) Joint modeling of multi-temporal stratospheric delay, deformation, and terrain residuals

[0144] Step 1: Build a multi-time relative spheric delay model window by window

[0145] Based on the window segmented in step (3), a linear model of tropospheric atmospheric delay phase and elevation is constructed pixel by pixel in any single window for a single interferogram:

[0146]

[0147] here It represents the tropospheric delay phase difference between the nth and n+1th SAR images and can be expressed as s is the total number of pixels in a single window, represents the tropospheric delay phase difference of the qth pixel (1≤q≤s); H hgt is the elevation parameter matrix corresponding to the pixels in the window, which can be expressed as h q Indicates the elevation value corresponding to the qth pixel (1≤q≤s); Representatives and The corresponding linear model parameter matrix of phase and elevation can be expressed as k n ,ξ n Represent the slope and intercept of the linear model of phase and elevation, respectively, and T represents the matrix or vector transpose.

[0148] Assume that InSAR data processing acquires N SAR images covering the same study area, and the acquisition time is (t1,…,t N ). A short baseline network is constructed according to the 36-day time baseline threshold and the 300m spatial baseline threshold, generating a total of M differential interferograms, M satisfying Here, it is assumed that the tropospheric delay of the reference image is zero. The tropospheric delay difference of two images at adjacent imaging moments can be used to represent the tropospheric delay of the M interferograms. Therefore, according to formula (1), the unwrapped phase of the M×s number of pixels in a single window of the M interferograms is obtained, which is specifically expressed as follows:

[0149]

[0150] here represents the unwrapped phase of M×s pixels in the same window of M interference patterns, where is the unwrapped phase of the qth pixel in the mth interferogram (1≤m≤M), The form is A represents the matrix related to the organization of the interference pattern, and A can be expressed as If the mth interferogram is obtained by interferometry of the i-th and j-th SAR images, then the corresponding element in the matrix A is A(m,l)=1 (i≤l<j), and the elements at other positions are zero; Indicates the calculation of the Kronecker product of two matrices; H hgt is the elevation parameter matrix corresponding to the s pixels in the window, which can be expressed as h q represents the elevation value corresponding to the qth pixel (1≤q≤s); X strat The linear model parameter matrix representing the tropospheric atmospheric delay phase and elevation of any adjacent SAR images can be expressed as k n ,ξ n represent the slope and intercept of the linear model of the tropospheric delay phase and elevation between the nth and n+1th SAR images (1≤n<N); They represent the deformation phase and terrain residual phase in the M interferograms, which will be used to construct the model in the second step; Δ represents the error vector, which mainly includes unwrapping error, residual turbulence, random noise, etc., and can be estimated through residual analysis or using an external error model.

[0151] Step 2: Construct surface deformation and terrain residual models window by window

[0152] Here, a cubic polynomial is used to describe the temporal variation of surface deformation. At the same time, the terrain residual is modeled based on the geometric model of the InSAR terrain phase. Therefore, the specific formulas for the deformation phase and terrain residual phase of each pixel in a single window can be expressed as follows:

[0153]

[0154] Here i , t j represent the imaging time of the i-th and j-th SAR images respectively; Respectively represent t i and t j The surface deformation phase and terrain residual phase of the qth pixel in time.

[0155] Δh q are the parameters to be estimated, representing the average velocity, acceleration, rate of change of acceleration, and terrain residual phase of the qth pixel respectively; λ is the radar wavelength, which can be taken as the 5.6cm wavelength of the Sentinel-1A radar satellite; r is the slant range between the SAR sensor and the ground, and θ is the incident angle of the SAR satellite sensor. Both r and θ can be obtained from the SLC parameter file; Indicates t i , t j The vertical baseline between time sensors can be calculated using the orbit parameter information of the Sentinel-1A satellite.

[0156] Then, formula (8) is extended to M interferograms, and the surface deformation phase and terrain residual phase of M×s pixels in a single window can be expressed as follows:

[0157]

[0158] Here I s×s represents the unit matrix of size s×s; Indicates the calculation of the Kronecker product of two matrices; the matrix G represents the time matrix corresponding to the matrix A, which can be expressed as M represents the total number of interferograms, N represents the total number of SAR images, and t n+1 -t n is the time interval between adjacent SAR images (1≤n<N). If the mth interferogram is obtained by interfering the i-th and j-th SAR images (i<j), then the corresponding element in the matrix G is G(m,n)=t n+1 -t n (1≤m≤M,n=i,…,j-1), the elements in other positions are zero;

[0159] D defo 、D topoRepresent the model coefficient matrices related to the surface deformation phase and the terrain residual phase, respectively. The specific form is expressed as Among them, t n+1 -t n is the interval between adjacent SAR images (1≤n<N), λ is the radar wavelength, which is taken here as the 5.6cm wavelength of the Sentinel-1A radar satellite. r is the slant range between the SAR sensor and the ground, and θ represents the SAR satellite sensor incidence angle. Both r and θ are obtained from the SLC parameter file. Indicates t i , t j The vertical baseline between time sensors can be calculated using the orbit parameter information of the Sentinel-1A satellite.

[0160] X defo 、X topo Represent the parameter vectors to be estimated, which can be expressed as Respectively represent the average velocity, acceleration, and acceleration change rate of the qth pixel; Δh q Represents the terrain residual phase of the qth pixel (1≤q≤s).

[0161] Step 3: Joint modeling of multi-temporal relativity spheric delay, surface deformation, and terrain residuals

[0162] In summary, by substituting formula (9) into formula (7), we can construct the observation equation of M interference patterns with M×s number of pixels in a single window:

[0163]

[0164] here is the actual observation value, representing the unwrapped phase of M×s pixels, expressed as B represents the coefficient matrix of the observation equation, expressed as X represents the parameter vector to be estimated, expressed as X=[X start ,X defo ,X topo ] T ; The meaning and representation of each symbol have been explained in the first and second steps; Δ represents the error vector of the observation equation, which mainly includes unwrapping error, residual turbulence, random noise, etc., and can be estimated through residual analysis or using an external error model.

[0165] (5) Estimation of model parameter M

[0166] The M-estimation is used to iteratively solve the parameter vector X to be estimated and the InSAR deformation time series in formula (10). The main steps include the following two parts:

[0167] Step 1: Calculate the initial weight matrix

[0168] The coherence γ of a single pixel g in the SAR interferogram m,g It can be obtained by calculating the mutual correlation coefficient of two complex SAR images, that is, calculating the overall mean of the complex signal sum of the local area centered on pixel g. The specific expression is:

[0169]

[0170] Here m,g represents the coherence of the g-th pixel in the m-th interference pattern (1≤m≤M,1≤g≤T), where T represents the total number of interference image elements; c i (μ) and c j (μ) represents the complex signal of each pixel in the local region centered on pixel g in the i-th and j-th SAR images, respectively, both of which are obtained from the SLC file. Here, i and j are the two SAR images used to generate the interferogram m (1≤i<j≤N); R is the number of pixels in the local region centered on pixel g; * is a conjugate complex number.

[0171] After calculating the coherence of each pixel in the interference pattern pixel by pixel using formula (11), the weight of the corresponding interference pattern can be determined by the following formula:

[0172]

[0173] Here m Indicates the weight of the mth interference pattern; represents the average coherence of the mth interference pattern and can be expressed as γ m,g is the coherence corresponding to each pixel calculated by formula (11) (1≤g≤T); C represents the average coherence of the unit weight interference pattern, which can be taken as 0.5 here.

[0174] Therefore, the initial weight matrix P can be expressed as follows:

[0175]

[0176] Step 2: M estimation algorithm flow

[0177] The solution of the parameter vector X under M estimation is obtained through formula (10). First, the ρ function is introduced as the loss function to optimize the parameters:

[0178]

[0179] ρ(x) is a real-valued function (not unique) and can be expressed as Here we usually take ω=1.345; represents the row vector, representing the ii-th row of the coefficient matrix B of the observation equation (1≤ii≤M×s); is the M-estimate solution of the parameter vector X; φ ii Represents the ii-th observation value, which is a column vector The unwrapped phase value of the value at the ii-th position in v ii is the observed value φ ii Correction number; p ii is the weight of the ii-th observation, which can be calculated from formula (13).

[0180] If the order After derivation of formula (14), the solution equation of M estimation is:

[0181]

[0182] Therefore, in the case of multiple parameters, the parameters The following iterative solution can be used:

[0183]

[0184] W is the equivalent weight matrix, which can be constructed using the ψ(x) function and the ρ(x) function, specifically expressed as W=diag(W1,W2,…,W M×s ),in or

[0185] In summary, M estimates the parameters The solution algorithm process is as follows:

[0186] ① Input data B, P, δ (where P is expressed as in formula (13); δ is a given threshold, which is 0.01 here);

[0187] ②Use the least squares method to find the initial value of the parameter

[0188] ③ Perform the following steps according to the process in ②:

[0189]

[0190] ④ Repeat steps ② and ③ until

[0191] ⑤Output results

[0192] (6) Phase correction and cumulative deformation time series solution

[0193] After completing the parameter estimation of the joint model of unwrapped phase and tropospheric delay, deformation and terrain residuals in a single window, the following steps will be continued:

[0194] ① Window-by-window parameter estimation: Execute steps (4) and (5) to complete the establishment and parameter estimation of the tropospheric delay model for all windows determined in step (3) one by one;

[0195] ② Parameter space interpolation and tropospheric delay reconstruction: First, the obtained tropospheric delay parameter values ​​are recorded at the positions of the geometric center pixels of each window. Then, the parameter values ​​at the geometric center positions of each window in each image sheet are used to interpolate the entire image sheet to obtain the parameter values ​​corresponding to all pixel positions in each image sheet. Finally, the tropospheric delay phase is reconstructed according to the elevation value corresponding to each pixel using formula (6). At this time, the tropospheric delay phase represents the tropospheric delay phase between SAR images at any adjacent imaging time.

[0196] ③ Convert the tropospheric delay phase: multiply the reconstructed tropospheric delay phase by a conversion matrix A (defined in formula (7)) to convert the tropospheric delay phase into a form corresponding to the interferogram connection mode;

[0197] ④ Unwrapped phase correction: remove the converted tropospheric delay phase and terrain residual phase from the unwrapped phase to complete the tropospheric atmospheric delay correction of the unwrapped phase;

[0198] ⑤ Deformation time series solution: Considering the relationship between vegetation density change and interferogram coherence, use coherence to quantify vegetation coverage, use the average coherence of the interferogram as the initial weight matrix, and apply the M estimation method in step (5) again to the unwrapped phase Φ′ after atmospheric delay correction. Replace the coefficient matrix B in formula (16) with matrix A, and the cumulative deformation phase within any adjacent SAR image imaging period can be iteratively solved.

[0199] Note: Generally, higher vegetation density indicates lower interferogram coherence during the corresponding period (existing literature confirms a negative correlation between vegetation cover (NDVI) and coherence), which can lead to more inaccurate or null values ​​in the interferogram. Here, coherence is used to quantify vegetation cover. By using the average coherence of each interferogram as an initial weighting factor in the calculation of the cumulative deformation phase, interferograms with lower coherence (high vegetation cover) are assigned lower weights to reduce their impact on the deformation results.

[0200] Figure 2 This is the system module connection diagram of the present invention. Figure 2As shown, the present invention also provides an InSAR atmospheric delay correction system based on adaptive window. Some modules of this system rely on GAMMA software. The system consists of seven modules: data preprocessing module, window adaptive dynamic segmentation module, InSAR data processing module, multi-time tropospheric delay, deformation and terrain residual joint modeling module, M estimation algorithm module, unwrapped phase tropospheric delay correction module and cumulative deformation time series solution module.

[0201] The data preprocessing module uses ERA5 reanalysis data to calculate the atmospheric water vapor content and its spatial location. Multi-temporal SAR data and DEM data preprocessing is mainly run on GAMMA software. The main preprocessing procedures include: SLC extraction, image registration, cropping, and multi-viewing.

[0202] The adaptive window dynamic segmentation module calculates the semivariogram function of atmospheric water vapor content and spatial location data using the variogram calculation module. A spherical or exponential model is then fitted to determine the model range. The number of SAR image segmentation windows is then determined based on the range. K-means clustering is used to classify longitude, latitude, and atmospheric water vapor content, resulting in multiple clusters of varying types. Finally, the SAR image window segmentation is completed by obtaining the longitude and latitude coordinates of each cluster boundary.

[0203] The InSAR data processing module calculates the interferometric baseline of the SLC product obtained by the data preprocessing module according to the set time baseline threshold and spatial baseline, and performs differential interferometry, filtering, phase unwrapping and elevation correction on the SAR image according to the connection of the baseline to obtain the unwrapped phase of the interferogram.

[0204] The multi-temporal relative stratosphere delay, deformation and terrain residual joint modeling module combines the segmentation window determined by the window adaptive dynamic segmentation module and the interferogram unwrapping phase obtained by the InSAR data processing module, and uses the elevation value, incident angle, slant range and vertical baseline parameters corresponding to each pixel in the window to construct a multi-temporal relative stratosphere delay, surface deformation and terrain residual joint model window by window.

[0205] The M estimation algorithm module offers two options for initial weight matrices: the first is the unit weight matrix. The second calculates the coherence between the pixels in the two SAR images corresponding to each interferogram based on baseline connectivity information from the InSAR data processing module and the complex signals of each pixel in the SAR image. Finally, the average coherence of each interferogram is calculated as the initial weight matrix. These initial weight matrices are selected based on different requirements, and the estimated parameters are solved through continuous iteration, achieving adaptive adjustment of residuals and ensuring the reliability of the estimation results.

[0206] The unwrapped phase tropospheric delay correction module first selects the unit weight matrix of the M estimation algorithm module as the initial weight matrix to iteratively solve the parameter model parameters, and then uses the model parameters to reconstruct the tropospheric delay phase and terrain residual phase. Finally, this part of the phase is removed from the unwrapped phase to complete the unwrapped phase tropospheric delay correction.

[0207] The cumulative deformation time series solution module uses the interferogram unwrapped phase corrected by the unwrapped phase tropospheric delay correction module as the observation value, and adopts the average coherence of the interferogram provided by the M estimation algorithm module as the initial weight matrix to iteratively solve the cumulative deformation time series.

[0208] By introducing an atmospheric moisture content product and adaptively segmenting SAR image windows based on differences in atmospheric water vapor content across space, this method effectively overcomes the loss of model accuracy caused by the traditional method's reliance on empirically set fixed window sizes. It also reduces computational redundancy during model parameter inversion. Furthermore, the relationship between vegetation cover density and interferogram coherence is considered, and the average coherence of the InSAR interferogram is used as the initial weight matrix to reduce the impact of seasonal vegetation density changes on deformation results. Furthermore, the M-estimation method is used to iteratively solve the InSAR deformation time series, improving the reliability of the deformation estimation results.

[0209] The above embodiments are not limitations of the present invention, and the present invention is not limited to the above examples. Any changes, modifications, additions or substitutions made by technicians in this technical field within the scope of the technical solution of the present invention also fall within the scope of protection of the present invention.

Claims

1. An InSAR atmospheric delay correction method based on an adaptive window, characterized by: The method comprises the following steps: Step S1, data acquisition and preprocessing; Step S2: InSAR data processing: Step S3, adaptive dynamic window segmentation: using the variogram to calculate the semivariogram function value list, the number of SAR image segmentation windows is determined by the size of the model range obtained by fitting, and the SAR image window segmentation is completed using the K-means clustering method; specifically, the following steps are included: Step 1: Calculate the variance function: Based on the spatial distribution of atmospheric moisture, use the variance function to calculate the semivariogram of atmospheric water vapor content and spatial position data, and select a spherical model or an exponential model for fitting to obtain the model range. Step 2: Window division and size adjustment: The number of SAR image segmentation windows is determined based on the range size. The longitude, latitude, and atmospheric water vapor content are classified using the K-means clustering method to obtain multiple clusters of different types. The longitude and latitude coordinates at the boundaries of each cluster are obtained to complete the SAR image window segmentation. Step S4, joint modeling of multi-temporal stratum delay, deformation and terrain residual, specifically includes the following steps: Step 1: Construct a multi-temporal relative stratigraphic delay model window by window; Step 2: Construct surface deformation and terrain residual models window by window; Step 3: Joint modeling of multi-temporal relativity spheric delay, surface deformation and topographic residuals; Step S5, estimating model parameters M; Step S6, unwrapping phase correction and cumulative deformation sequence solution: the average coherence of the interferogram is introduced as the initial weight matrix, and the M estimation method is used again to solve the deformation time series to achieve atmospheric delay correction of the deformation results.

2. The InSAR atmospheric delay correction method based on adaptive window according to claim 1, characterized in that: In step S1, the following data need to be acquired: multi-temporal InSAR data, ERA5 data corresponding to the time when the InSAR data was acquired, and DEM data; the above data are preprocessed: Multi-temporal InSAR data preprocessing includes SLC extraction, DEM geocoding, image registration, study area cropping, and multi-viewing; ERA5 data preprocessing includes data format conversion, spatial region cropping, vertical integration of temperature and relative humidity to obtain the atmospheric water vapor content of the entire troposphere, and then calculating the average atmospheric water vapor content in the spatial region and performing spatial interpolation for insufficient spatial resolution.

3. The InSAR atmospheric delay correction method based on adaptive window according to claim 2, characterized in that: In step S2, the pre-processed SLC image is first subjected to interferometric baseline calculation according to the time baseline threshold and the spatial baseline threshold. Then, based on the connection of the baselines, the SAR image is subjected to differential interferometry, filtering, phase unwrapping, and elevation correction. The phase formula after unwrapping is as follows: （2）， in, represents the unwrapped phase of the entire interferogram, Represents the deformation phase, represents the terrain residual phase, represents the tropospheric delayed phase, represents the error vector.

4. The InSAR atmospheric delay correction method based on adaptive window according to claim 3, characterized in that: The step 1 specifically includes the following steps: Based on the average atmospheric water vapor content and spatial location data obtained from ERA5 data, the semivariogram function is used as the variation function to describe the degree of change in water vapor content between different locations. The specific calculation formula is: （2）， in, is the semivariogram function value, and In spatial position 、 The atmospheric water vapor content at is the spatial distance, The spatial distance is The total number of pixel pairs; Use the Euclidean distance formula to calculate the distance between each pixel pair (x i ,y i ) and (x j ,y j ) Calculate spatial distance : （3）， According to the calculated spatial distance , classify all pixel pairs, group the pixel pairs within the same distance range into one category, and for each category, the same Calculate the semivariogram function value according to formula (2) and record the distance between each pixel. The corresponding semivariogram function values ​​form a list of semivariogram function values.

5. The InSAR atmospheric delay correction method based on adaptive window according to claim 4, characterized in that: The step 2 specifically includes the following steps: Use the semivariogram value list calculated in step 1 to determine the window size, select an appropriate semivariogram model for fitting, and obtain the model range. According to the range value and the range of the SAR image processed in step S2, that is, the size of the SLC image obtained by data preprocessing, the number of windows in the spatial area covered by the SLC image is preliminarily determined. Then, the longitude, latitude, and water vapor content values ​​of each pixel in the SAR image are used as the sample set, and the sample set is divided into multiple clusters according to the distance between samples. , each cluster must satisfy the square error The value of is the smallest: ， in, Indicates the Clusters The sample value vector corresponding to each pixel in ; Indicates the Clusters The mean vector of , called the centroid, is expressed as , It is The number of samples in a cluster; After the sample set is divided, the centroid and longitude and latitude values ​​at the boundary of each cluster are extracted, and the coordinate vector of the inflection point at the boundary of each window is determined according to the following formula: （5）， in, To ensure the The new boundary inflection point coordinate vector is calculated based on the 20% overlap area between the cluster and the adjacent cluster; It means that the first The coordinate vector of the inflection point at the cluster boundary; Indicates the The coordinates of the centroid position of each cluster; In summary, the coordinate vector of the boundary inflection point of each cluster The SAR image segmentation can be completed.

6. The InSAR atmospheric delay correction method based on adaptive window according to claim 3, characterized in that: The step S4 specifically includes the following steps: Step 1: Based on the window segmented in step S3, a linear model of tropospheric atmospheric delay phase and elevation is constructed pixel by pixel in any single window for a single interferogram: （6）， in, Indicates the Frame and The tropospheric delay phase difference between the SAR images is is the elevation parameter matrix corresponding to the pixels in the window, Representatives and The corresponding linear model parameter matrix of phase and elevation; is the total number of pixels in a single window, Indicates the The tropospheric delay phase difference of pixels, ; Indicates the The elevation value corresponding to the pixel, ; 、 represent the slope and intercept of the linear model of phase and elevation, respectively; ， ， ， InSAR data processing has obtained a total of SAR images, generating Amplitude differential interferogram, satisfy , the tropospheric delay of the reference image is zero, according to formula (1) we get Interference patterns in a single window The unwrapped phase of the number of pixels is expressed as follows: （7）， in, express Interference patterns in the same window The unwrapped phase of the number of pixels, represents the matrix related to the organization of the interference pattern, Indicates the calculation of the Kronecker product of two matrices. For this window The elevation parameter matrix corresponding to the pixels, The linear model parameter matrix representing the tropospheric atmospheric delay phase and elevation of any adjacent SAR images, Respectively The deformation phase and topographic residual phase in the interferogram will be used to construct the model in the second step. represents the error vector, estimated by residual analysis or using an external error model; For the Interference pattern The unwrapped phase of pixels, , Indicates the The elevation value corresponding to the pixel, 、 Representing the Frame and The slope and intercept of the linear model of the tropospheric delay phase and elevation between SAR images, ; ， ， ， ； Step 2: Use a cubic polynomial to describe the temporal variation of surface deformation. Simultaneously, model the terrain residual based on the geometric model of the InSAR terrain phase. The specific formulas for the pixel-by-pixel deformation phase and terrain residual phase within a single window are as follows: （8）， in, 、 Representing the 、 Imaging time of a SAR image; Respectively and First The surface deformation phase and terrain residual phase of each pixel; are the parameters to be estimated, representing the The average velocity, acceleration, acceleration change rate and terrain residual phase of each pixel; is the radar wavelength; Indicates the slant distance between the SAR sensor and the ground, is the SAR satellite sensor incident angle, 、 All are obtained from the SLC parameter file; express 、 vertical baseline between time sensors; Expand formula (8) to Interference pattern, then within a single window The surface deformation phase and terrain residual phase of the number pixel are expressed as follows: （9）， in, express The unit matrix of size; matrix Representatives and matrices The corresponding time matrix is ​​expressed as , represents the total number of interference patterns, represents the total number of SAR images, is the time interval between adjacent SAR images, ; 、 They represent the model coefficient matrices related to the surface deformation phase and the terrain residual phase, respectively. The specific form is expressed as: ， in, Representative The terrain residual phase of pixels, ; 、 Represent the parameter vectors to be estimated respectively; ； Step 3: Substitute formula (9) into formula (7) to construct Interference patterns in a single window The observation equation of the number of pixels: （10）， in, is the actual observed value, representing The unwrapped phase of the number of pixels, represents the coefficient matrix of the observation equation, Represents the parameter vector to be estimated, expressed as ; 。 7. The InSAR atmospheric delay correction method based on adaptive window according to claim 6, characterized in that: In step S5, the parameter vector to be estimated in formula (10) is estimated using M. The iterative solution of the InSAR deformation time series includes the following two steps: Step 1: Calculate the initial weight matrix Calculate pixel The overall mean of the complex signal sum of the local area centered is expressed as: （11）， in, Indicates the Interference pattern The coherence of pixels, , represents the total number of interference image elements; and Respectively represent Frame and Pixels in a SAR image The complex signals of each pixel in the local area centered are obtained from the SLC file. and Is used to generate the interference pattern Two SAR images ; Pixel The number of pixels in the local area centered on ; is a conjugate complex number; After calculating the coherence of each pixel in the interference pattern pixel by pixel using formula (11), the weight of the corresponding interference pattern is determined by the following formula: （12）， in, Indicates the The weight of the interference pattern; Representative The average coherence of the interference pattern is expressed as , represents the average coherence of the unit weight interference pattern; Therefore, the initial weight matrix P is expressed as: （13）， Step 2: M estimation algorithm process The parameter vector is obtained by formula (10) The solution under M estimation is first introduced The function is used as the loss function to optimize the parameters: （14）， in, is a real-valued function and is not unique, expressed as , ; Represents a row vector, representing the coefficient matrix of the observation equation No. OK, ; is the parameter vector The M-estimate solution of ; is the observed value number of corrections; This is for the The weight of each observation is calculated from formula (13); Indicates the observations, a column vector The The unwrapped phase value of the value at each position; The column vector is represented as ; make , after derivation of formula (14), the solution equation of M estimation is: （15）， Therefore, in the case of multiple parameters, the parameters Solve it iteratively using the following formula: （16）， in, is the equivalent weight matrix, which uses Function and Function construction, specifically expressed as ,in .

8. The InSAR atmospheric delay correction method based on adaptive window according to claim 7, characterized in that: In step S6, after completing the parameter estimation of the joint model of the unwrapped phase and tropospheric delay, deformation and terrain residual in a single window, the following steps will be continued: S61, window-by-window parameter estimation: executing steps S4 and S5 to complete the establishment and parameter estimation of the tropospheric delay model for all windows determined in step S3 one by one; S62. Parameter space interpolation and tropospheric delay reconstruction: First, the obtained tropospheric delay parameter values ​​are recorded at the positions of the geometric center pixels of their respective windows. Then, the parameter values ​​at the geometric center positions of each window in each image sheet are used to interpolate the entire image sheet to obtain the parameter values ​​corresponding to all pixel positions in each image sheet. Finally, the tropospheric delay phase is reconstructed according to the elevation value corresponding to each pixel using formula (6). S63, transform the tropospheric delay phase: multiply the reconstructed tropospheric delay phase by a transformation matrix To convert the tropospheric delayed phase into a form corresponding to the interferogram connection mode; S64, unwrapped phase correction: removing the converted tropospheric delay phase and terrain residual phase from the unwrapped phase to complete the tropospheric atmospheric delay correction of the unwrapped phase; S65, Deformation time series solution: Use the average coherence of the interferogram as the initial weight matrix to correct the unwrapped phase after atmospheric delay Apply the M estimation method of step S5 again: transform the coefficient matrix in formula (16) into Replace with matrix , iteratively solve the cumulative deformation phase in any adjacent SAR image imaging period .

9. A system for implementing the adaptive window-based InSAR atmospheric delay correction method according to any one of claims 1 to 8, characterized in that: The system includes: Data preprocessing module: uses ERA5 data to calculate the atmospheric water vapor content and its spatial location; Window adaptive dynamic segmentation module: SAR image window segmentation is completed based on the value list obtained by the data preprocessing module according to the variance function calculation and combined with the K-means clustering method; InSAR data processing module: The interferometric baseline is calculated for the SLC image obtained by the data preprocessing module according to the set time baseline threshold and spatial baseline. Based on the connection of the baseline, the SAR image is subjected to differential interferometry, filtering, phase unwrapping, and elevation correction to obtain the interferogram unwrapping phase. Multi-temporal relative spheric delay, deformation, and terrain residual joint modeling module: This module combines the segmentation window determined by the window adaptive dynamic segmentation module and the interferogram unwrapped phase obtained by the InSAR data processing module, and uses the elevation value, incident angle, slant range, and vertical baseline corresponding to each pixel in the window to construct a multi-temporal relative spheric delay, surface deformation, and terrain residual joint model window by window; M estimation algorithm module: This module provides two initial weight matrices. The first is the unit weight matrix. The second is to calculate the coherence between the corresponding pixels of the two SAR images corresponding to each interferogram based on the baseline connection information in the InSAR data processing module and the complex signal of each pixel in the SAR image. Finally, the average coherence of each interferogram is calculated as the initial weight matrix. The estimated parameters are solved through continuous iteration to achieve adaptive adjustment of the residual. Unwrapped phase tropospheric delay correction module: This module first uses the unit weight matrix of the M estimation algorithm module as the initial weight matrix to iteratively solve the model parameters, then uses the model parameters to reconstruct the tropospheric delay phase and terrain residual phase, and finally removes this part of the phase from the unwrapped phase to complete the unwrapped phase tropospheric delay correction; Cumulative deformation time series solution module: The unwrapped phase of the interferogram corrected by the unwrapped phase tropospheric delay correction module is used as the observation value, and the average coherence of the interferogram provided by the M estimation algorithm module is used as the initial weight matrix to iteratively solve the cumulative deformation time series.

Images