Method, system and computer readable medium for atmospheric delay correction of time series insar monitoring data

Abstract

The application provides a kind of atmospheric delay correction method, system and computer readable medium of time series InSAR monitoring data, which can serve the accurate estimation and correction of InSAR atmospheric delay in high mountain and canyon area, realize the reading of time series InSAR high coherence point correlation information based on the input SAR image, such as phase, elevation, initial deformation rate, unwrapping closed loop error, etc., window segmentation is carried out on InSAR interferogram by sliding multi-window with a certain overlap degree, linear model parameters in each window are iteratively estimated, and finally the estimated average deformation rate and time series deformation variable of each point are output. Through the application, the atmospheric delay correction of InSAR monitoring data in high mountain and canyon area dominated by vertical stratification can be realized, which provides an effective way for time series InSAR technology in high mountain and canyon area to overcome the influence of atmospheric delay, and improves the precision and applicability of time series InSAR technology in surface deformation monitoring in high mountain and canyon area.

Description

Technical Field

[0001] This invention relates to the field of InSAR monitoring technology, and more specifically to an atmospheric delay correction method, system, and computer-readable medium for time-series InSAR monitoring data. Background Technology

[0002] Frequent landslides and other geological disasters in mountainous areas cause severe economic and human losses. Identifying, monitoring, and providing early warning of potential landslide hazards in high mountain and canyon areas poses a significant challenge to regional development. Temporal InSAR technology, with its advantages of all-day, all-weather, and wide-area coverage, has been applied to deformation monitoring in landslides, earthquakes, mining areas, and cities, generating substantial social and economic benefits. However, in high mountain and canyon areas, temporal InSAR is severely affected by the seasonal vertical stratification delay within the tropospheric delay, significantly interfering with deformation monitoring results. This delay is highly correlated with topography. Summary of the Invention

[0003] In view of the defects and shortcomings of the existing technology, the first aspect of the present invention proposes an atmospheric delay correction method for time-series InSAR monitoring data, comprising:

[0004] (1) After fusing distributed scatterers and permanent scatterers using DS-InSAR technology and performing phase stability analysis, the entangled phase of highly coherent measurement points is obtained.

[0005] (2) Perform 3D phase unwrapping to obtain the unwrapped phase values ​​of the measurement points of each interferogram;

[0006] (3) Based on the preset window size and window movement step, each interferogram is divided into independent windows to obtain the point information, unwrapping phase, initial deformation rate and elevation value in each window;

[0007] (4) Within a certain window i, stable points are selected based on the initial deformation rate or the updated average deformation rate. Then, a linear model equation of phase and elevation is established on the selected stable points. The least squares method is used to solve the linear model parameters of each phase in the window. The tropospheric delay phase of all measurement points in the window is calculated based on the elevation value.

[0008] (5) After subtracting the tropospheric delay phase calculated for the first time from the initial unwrapped phase, the atmospheric delay correction phase is obtained. The deformation rate of each point is re-estimated using the weighted least squares method. The atmospheric delay phase is estimated again by repeating the process of step (4).

[0009] (6) Calculate the root mean square of the difference between two adjacent estimates of atmospheric delay phase for each interferogram within window i, and calculate the average value of all interferograms; determine the termination threshold after at least three estimates of atmospheric delay phase have been completed. Is it less than a preset threshold? If yes, proceed to the next window estimation; if no, return to step (4) to continue estimating the atmospheric delay phase.

[0010] (7) After the atmospheric delay phase is estimated in all windows, the biharmonic spline interpolation method is applied to interpolate the linear model parameters at the center of each window to each measurement point, and the atmospheric phase delay of each measurement point is calculated in combination with the elevation of each point.

[0011] (8) Subtract the atmospheric phase delay of the measurement point from the winding phase of the measurement point, perform 3D phase unwrapping again, and further determine whether the change value ΔU of the percentage of error-free points in the two unwrappings is less than the threshold setting. If so, stop the iteration; otherwise, return to step (3).

[0012] (9) After atmospheric delay correction is completed, the residual turbulence delay is removed by spatiotemporal filtering method;

[0013] (10) The average deformation rate and time series deformation of each point after atmospheric correction are obtained by modeling deformation information.

[0014] According to a second aspect of the present invention, a computer system is also provided, comprising:

[0015] One or more processors;

[0016] A memory that stores operable instructions;

[0017] The instructions, when executed by the one or more processors, implement the aforementioned atmospheric delay correction method for time-series InSAR monitoring data.

[0018] According to a third aspect of the invention, a computer-readable medium storing a computer program is also provided, the computer program including instructions executable by one or more computers, the instructions, when executed by the one or more computers, performing the aforementioned atmospheric delay correction method for time-series InSAR monitoring data.

[0019] As can be seen from the above technical solutions of the present invention, the atmospheric delay correction method for time-series InSAR monitoring data proposed in this invention aims to address the problem that time-series InSAR is severely affected by the seasonal vertical stratification delay in the tropospheric delay, which seriously interferes with deformation monitoring results. This invention proposes an accurate estimation and correction method for atmospheric delay in InSAR in high mountain and canyon areas. Specifically, the atmospheric delay correction method for time-series InSAR monitoring data obtains the temporal winding phase of high-density measurement points based on time-series InSAR technology and proposes an optimized multi-temporal move-window linear model (OMMLM) suitable for high mountain and canyon areas. This achieves atmospheric delay correction for InSAR monitoring data in high mountain and canyon areas dominated by vertical stratification, providing an effective way for time-series InSAR technology in high mountain and canyon areas to overcome the influence of atmospheric delay and improve the accuracy and applicability of time-series InSAR technology in monitoring surface deformation in high mountain and canyon areas.

[0020] It should be understood that all combinations of the foregoing concepts and the additional concepts described in more detail below may be considered part of the inventive subject matter of this disclosure, provided that such concepts do not contradict each other. Furthermore, all combinations of the claimed subject matter are considered part of the inventive subject matter of this disclosure.

[0021] The foregoing and other aspects, embodiments, and features of the teachings of the present invention will be more fully understood from the following description in conjunction with the accompanying drawings. Other additional aspects of the invention, such as features and / or beneficial effects of exemplary embodiments, will become apparent from the following description or may be learned through practice of specific embodiments according to the teachings of the present invention. Attached Figure Description

[0022] The accompanying drawings are not intended to be drawn to scale. In the drawings, each identical or nearly identical component shown in the various figures may be denoted by the same reference numeral. For clarity, not every component is labeled in each figure. Embodiments of various aspects of the invention will now be described by way of example and with reference to the accompanying drawings.

[0023] Figure 1 This is a schematic diagram illustrating the principle of the atmospheric delay correction method for time-series InSAR monitoring data according to an embodiment of the present invention.

[0024] Figure 2 This is a flowchart illustrating the atmospheric delay correction method for time-series InSAR monitoring data according to an embodiment of the present invention.

[0025] Figure 3 This is a schematic diagram of the movement of a rectangular window according to an embodiment of the present invention.

[0026] Figure 4This is an unwrapped interferogram after atmospheric delay correction according to an embodiment of the present invention. Detailed Implementation

[0027] To better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0028] Various aspects of the invention are described in this disclosure with reference to the accompanying drawings, which illustrate numerous illustrative embodiments. The embodiments of this disclosure are not necessarily intended to encompass all aspects of the invention. It should be understood that the various concepts and embodiments described above, as well as those described in more detail below, can be implemented in any of many ways, because the concepts and embodiments disclosed herein are not limited to any particular implementation. Furthermore, some aspects of the invention disclosed may be used alone or in any suitable combination with other aspects of the invention disclosed.

[0029] The atmospheric delay correction method for time-series InSAR monitoring data, combined with the embodiments of the present invention, can be mainly used for accurate estimation and correction of InSAR atmospheric delay in high mountain and canyon areas. It can read relevant information of high coherence points in time-series InSAR based on the input SAR image, such as phase, elevation, initial deformation rate, unwrapping and closure error, etc. The InSAR interferogram is divided into windows by a sliding multi-window with a certain degree of overlap, and the linear model parameters in each window are estimated iteratively. Finally, the estimated average deformation rate and time series deformation of each point are output.

[0030] Combination Figure 1 , 2 As shown, the atmospheric delay correction method for time-series InSAR monitoring data, as an example, includes the following steps:

[0031] (1) After fusing distributed scatterers and permanent scatterers using DS-InSAR technology and performing phase stability analysis, the entangled phase of highly coherent measurement points is obtained.

[0032] (2) Perform 3D phase unwrapping to obtain the unwrapped phase values ​​of the measurement points of each interferogram;

[0033] (3) Based on the preset window size and window movement step, each interferogram is divided into independent windows to obtain the point information, unwrapping phase, initial deformation rate and elevation value in each window;

[0034] (4) Within a certain window i, stable points are selected based on the initial deformation rate or the updated average deformation rate. Then, a linear model equation of phase and elevation is established on the selected stable points. The least squares method is used to solve the linear model parameters of each phase in the window. The tropospheric delay phase of all measurement points in the window is calculated based on the elevation value.

[0035] (5) After subtracting the tropospheric delay phase calculated for the first time from the initial unwrapped phase, the atmospheric delay correction phase is obtained. The deformation rate of each point is re-estimated using the weighted least squares method. The atmospheric delay phase is estimated again by repeating the process of step (4).

[0036] (6) Calculate the root mean square of the difference between two adjacent estimates of atmospheric delay phase for each interferogram within window i, and calculate the average value of all interferograms; determine the termination threshold after at least three estimates of atmospheric delay phase have been completed. Is it less than a preset threshold? If yes, proceed to the next window estimation; if no, return to step (4) to continue estimating the atmospheric delay phase.

[0037] (7) After the atmospheric delay phase is estimated in all windows, the biharmonic spline interpolation method is applied to interpolate the linear model parameters at the center of each window to each measurement point, and the atmospheric phase delay of each measurement point is calculated in combination with the elevation of each point.

[0038] (8) Subtract the atmospheric phase delay of the measurement point from the winding phase of the measurement point, perform 3D phase unwrapping again, and further determine whether the change value ΔU of the percentage of error-free points in the two unwrappings is less than the threshold setting. If so, stop the iteration; otherwise, return to step (3).

[0039] (9) After atmospheric delay correction is completed, the residual turbulence delay is removed by spatiotemporal filtering method;

[0040] (10) The average deformation rate and time series deformation of each point after atmospheric correction are obtained by modeling deformation information.

[0041] In the embodiments of this invention, it should be understood that Satellite Synthetic Aperture Radar Interferometry (InSAR) is a radar technology applicable to mapping and remote sensing. It uses a satellite or aircraft-borne synthetic aperture radar system to acquire high-resolution complex images of ground reflections. Each resolution pixel contains not only grayscale information but also the phase signal required for interferometry. By using two or more synthetic aperture radar images, digital elevation models or surface deformation maps can be generated based on the phase difference of the echoes received by the satellite or aircraft. This can also be used for monitoring natural disasters such as earthquakes, volcanoes, and landslides, as well as structural engineering, especially settlement monitoring and structural stability.

[0042] Distributed Target Radar Interferometry (DS-InSAR) is developed based on InSAR technology. It uses SAR data acquired at different times in the same area for differential interferometric processing. By removing common factors (flatland effect, terrain phase, atmospheric delay, etc.) from the two observation phases through differential processing, the deformation phase is obtained, thereby obtaining surface deformation information. DInSAR technology inherits the advantages of InSAR technology, such as wide observation range and high degree of automation. It can also obtain more temporal InSAR coherence point density in complex and dangerous mountainous environments. In landslide monitoring, it can estimate the subsidence boundary range, analyze the subsidence evolution process, and obtain the overall activity range and intensity of landslides, which has significant advantages.

[0043] In Distributed Target Radar Interferometry (DS-InSAR) technology, a scatterer refers to an object within a pixel that may contain multiple reflected signals. These are mainly categorized into two types: distributed scatterers and permanent scatterers (also known as persistent scatterers, PS). Permanent scatterers typically refer to pixels containing objects whose reflectivity is time-coherent, remaining constant throughout the entire time period. Examples include fixed buildings such as houses in urban areas, and man-made structures such as iron towers.

[0044] Distributed scatterers (DS) typically refer to pixels in which multiple objects exhibit weak reflections with relatively similar intensities. Generally, a distributed scatterer indicates that the individual scattered echo intensities of multiple objects within a single resolution cell are not very high, but within a certain range, they follow the same statistical distribution as their surrounding adjacent targets. Their echo signal scattering model statistically satisfies a Gaussian distribution, manifesting as the vector sum of multiple echoes. Distributed scatterers generally correspond to land cover types with low vegetation cover, such as bare land and deserts.

[0045] In an embodiment of the present invention, phase stability analysis is performed based on the existing DS-InSAR technology, which integrates the settings of distributed scatterers and permanent scatterers, to obtain the entangled phase of highly coherent measurement points.

[0046] As an optional example, in the embodiments of the present invention, the SAR image data reading and analysis are performed using the MATLAB-based InSAR software Stamps (v4.1 version as an example). High coherence (high density) monitoring points can be obtained, and the unwrapped phase, elevation value change rate and unwrapped closed-loop error parameters can be obtained through phase unwrapping.

[0047] For ease of explanation, assume that M pairs of interferograms are generated from N SAR images, and assume that the tropospheric delay of the reference super master image selected during registration is zero, and that residual terrain and orbital errors have been removed from each interferogram.

[0048] As an optional embodiment, in step (3), each interferogram is independently segmented using a rectangular moving window with overlap, wherein the width of the rectangular moving window in the azimuth / range direction and the window sliding step size are respectively represented by S. windows and S steps And S windows =2S steps .

[0049] Combination Figure 2 , 3 As shown, each interferogram is divided into independent windows by sliding the window, resulting in I windows. The number of windows I is calculated using the following formula:

[0050]

[0051] Among them, S range S represents the number of range pixels in the interferogram; azimuth For the number of azimuth pixels in the interferogram, This represents the floor operation.

[0052] According to the phase-elevation linear model (LM), within a small local area of ​​each window, the tropospheric delay phase and elevation dominated by vertical stratification satisfy the following linear relationship:

[0053]

[0054] Where, φ i,n This represents the tropospheric delay phase at any point within the i-th moving window of the n-th SAR image acquisition time; n = 1, 2, ..., N, i = 1, 2, ..., I, k i,n The slope reflects the scale of the delay in the vertical stratification of the troposphere; h i,n This represents the elevation value, in kilometers. The intercept represents the large-scale tropospheric component and a portion of the turbulent delay component.

[0055] To reduce the impact of seasonal variations in local stochastic turbulence mixing delay and vertical stratification delay on parameter estimation, this method utilizes multi-temporal observation phases to ensure the robustness of local window linear estimation. The unwrapped phase of J points within the i-th moving window of M interferograms is specifically expressed as follows:

[0056]

[0057] Among them, H i Represents the elevation matrix; A i The linear transformation coefficient matrix represents the transformation from SAR image to interferometric image pair, recording the main and auxiliary image combination for each pair of interferograms; X represents the tensor product;i Represents the matrix of unknown parameters to be solved in the linear model, including the slope and intercept of the linear model at each time phase; E is a unit column vector; Φ i RES is a column vector consisting of M interferometric unwrapping phase values ​​at each point within the window; i This represents the residual phase after removing terrain and orbital errors, including deformation signals, residual turbulence, and random noise.

[0058] To reduce the impact of the deformation region on parameter estimation, in step (4), points whose absolute value of deformation rate is less than an empirical threshold are selected as stable points and weights are assigned. A weight matrix P is established based on the initial deformation rate, expressed in the following form:

[0059] P = 1 / (1 + |v) sta |),|v sta |≤v thr

[0060] Among them, v sta The deformation rate at the steady point; v thr The empirical threshold for deformation rate is set to 5 mm / a.

[0061] Combination Figure 2 As can be seen, in the correction method of the embodiment of the present invention, after the first estimation of the tropospheric delay phase within the window range, the preliminary tropospheric delay phase of each point is obtained. The obtained tropospheric delay phase is subtracted from the unwrapped phase and the deformation rate of each point is re-estimated. Stable points are re-selected based on the new deformation rate until the difference between the tropospheric delay phases estimated by two adjacent times is within the expected range and there is no significant difference.

[0062] As an optional embodiment, in step (6), a termination threshold is calculated for the atmospheric delay phase difference estimated before and after the same window. Furthermore, it determines the relationship between the value and a preset threshold. If the value is less than the preset threshold, the next window is estimated; otherwise, it returns to step (4) to continue estimating the atmospheric delay phase of the current window.

[0063] First, define the root mean square of the two atmospheric delay phase parameter estimates for a given window.

[0064]

[0065] in, This represents the difference in tropospheric delay phase between two estimates of the m-th interferometric pair point j; M is the number of interferograms; J is the number of measurement points within the window;

[0066] Then, to Perform a difference operation to obtain the termination threshold of the new window loop.

[0067]

[0068] That is, at least three window parameter estimations must be completed, when two consecutive windows are in use. When the difference is less than a preset threshold, it indicates that the degree of change in the difference between two adjacent parameter estimation results meets expectations, and further increasing the number of cycles can no longer significantly improve the atmospheric delay estimation results, meaning that the tropospheric delay model parameter estimation results have tended to stabilize at this point. Figure 4 As shown, an example is an unwrapped interferogram after atmospheric correction.

[0069] As an optional embodiment, in step (8), after subtracting the atmospheric delay phase from the initial wrapped phase before atmospheric delay correction, 3D phase unwrapping is performed again, and the change in the percentage of error-free points between the two unwrapping operations relative to the total number of interferogram measurement points is calculated as ΔU:

[0070]

[0071] Among them, U i+1 and U i N represents the percentage of error-free unwrapped points in all interferograms after the (i+1)th and ithth estimations of the atmospheric delay phase and the removal of the interference phase, respectively; mp This indicates the number of measurement points; M is the number of interferograms; error-free points refer to points where the absolute value of the unwrapped closed-loop residual is less than 1 rad.

[0072] According to an embodiment of the present invention, a computer system is also provided, comprising: one or more processors and a memory storing operable instructions. The aforementioned instructions, when executed by the one or more processors, are used to implement the aforementioned atmospheric delay correction method for time-series InSAR monitoring data.

[0073] According to an embodiment of the present invention, a computer-readable medium for storing a computer program is also provided, the computer program including instructions executable by one or more computers, which, when executed by one or more computers, complete the process of the aforementioned atmospheric delay correction method for time-series InSAR monitoring data.

[0074] While the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. An atmospheric delay correction method for time-series InSAR monitoring data, characterized in that, include: (1) After fusing distributed scatterers and permanent scatterers using DS-InSAR technology and performing phase stability analysis, the entangled phase of highly coherent measurement points is obtained. (2) Perform 3D phase unwrapping to obtain the unwrapped phase values ​​of the measurement points of each interferogram; (3) Based on the preset window size and window movement step, each interferogram is divided into independent windows to obtain the point information, unwrapping phase, initial deformation rate and elevation value in each window; (4) Within a certain window i, stable points are selected based on the initial deformation rate or the updated average deformation rate. Then, a linear model equation of phase and elevation is established on the selected stable points. The least squares method is used to solve the linear model parameters of each phase in the window. The tropospheric delay phase of all measurement points in the window is calculated based on the elevation value. (5) After subtracting the tropospheric delay phase calculated for the first time from the initial unwrapped phase, the atmospheric delay correction phase is obtained. The deformation rate of each point is re-estimated using the weighted least squares method. The atmospheric delay phase is estimated again by repeating the process of step (4). (6) Calculate the root mean square of the difference between two adjacent estimates of atmospheric delay phase for each interferogram within window i, and calculate the average value of all interferograms; determine the termination threshold after at least three estimates of atmospheric delay phase have been completed. Is it less than a preset threshold? If yes, proceed to the next window of estimation; if no, return to step (4) to continue estimating the atmospheric delay phase; where the above-mentioned termination threshold is... Indicates two consecutive times difference, This represents the root mean square of the estimates of the atmospheric delay phase parameters for a given window. (7) After the atmospheric delay phase is estimated in all windows, the biharmonic spline interpolation method is applied to interpolate the linear model parameters at the center of each window to each measurement point, and the atmospheric phase delay of each measurement point is calculated in combination with the elevation of each point. (8) Subtract the atmospheric phase delay of the measurement point from the winding phase of the measurement point, perform 3D phase unwrapping again, and further determine whether the change value ΔU of the percentage of error-free points in the two unwrappings is less than the threshold setting. If so, stop the iteration; otherwise, return to step (3). (9) After atmospheric delay correction is completed, the residual turbulence delay is removed by spatiotemporal filtering method; (10) The average deformation rate and time series deformation of each point after atmospheric correction are obtained by modeling deformation information.

2. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 1, characterized in that, In step (3), each interferogram is independently segmented using a rectangular moving window with overlap, wherein the width of the rectangular moving window in the azimuth / range direction and the window sliding step size are respectively represented by S. windows and S steps And S windows =2S steps .

3. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 2, characterized in that, In step (3), I windows are obtained by dividing each interferogram into independent windows. The number of windows is calculated using the following formula: Among them, S range S represents the number of range pixels in the interferogram; azimut h represents the number of azimuth pixels in the interferogram. This represents the floor operation; Within each window, the tropospheric delay phase and elevation dominated by vertical stratification satisfy the following linear relationship: Where, φ i,n This represents the tropospheric delay phase at any point within the i-th moving window of the n-th SAR image acquisition time; n = 1, 2, ..., N, i = 1, 2, ..., I, k i,n The slope reflects the scale of the delay in the vertical stratification of the troposphere; h i,n This represents the elevation value, in kilometers. The intercept represents the large-scale tropospheric component and a portion of the turbulent delay component.

4. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 3, characterized in that, The 3D phase unwrapping includes: The unwrapping phase expression for J points within the i-th moving window of M interferograms is as follows: Among them, H i Represents the elevation matrix; A i The linear transformation coefficient matrix represents the transformation from SAR image to interferometric image pair, recording the main and auxiliary image combination for each pair of interferograms; X represents the tensor product; i Represents the matrix of unknown parameters to be solved in the linear model, including the slope and intercept of the linear model at each time phase; E is a unit column vector; Φ i RES is a column vector consisting of M interferometric unwrapping phase values ​​at each point within the window; i This represents the residual phase after removing terrain and orbital errors, including deformation signals, residual turbulence, and random noise.

5. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 1, characterized in that, In step (4), points whose absolute value of deformation rate is less than an empirical threshold are selected as stable points and weights are assigned. A weight matrix P is established based on the initial deformation rate, expressed in the following form: P=1 / (1+|v sta |),|v sta |≤v thr ; Among them, v sta The deformation rate at the steady point; v thr The empirical threshold for deformation rate is set to 5 mm / a.

6. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 1, characterized in that, In step (5), after the first tropospheric delay phase estimation within the window range, the preliminary tropospheric delay phase of each point is obtained. The obtained tropospheric delay phase is subtracted from the unwrapped phase and the deformation rate of each point is re-estimated. Stable points are re-selected based on the new deformation rate until the difference between the two adjacent estimated tropospheric delay phases is within the expected range.

7. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 1, characterized in that, In step (6), a termination threshold is calculated for the atmospheric delay phase difference estimated before and after the same window. Furthermore, it determines the relationship between the value and a preset threshold. If the value is less than the preset threshold, the next window is estimated; otherwise, it returns to step (4) to continue estimating the atmospheric delay phase of the current window. First, define the root mean square of the two atmospheric delay phase parameter estimates for a given window. in, This represents the difference in tropospheric delay phase estimated before and after point j of the m-th interferometric pair; M is the number of interferograms; J is the number of measurement points within the window; Then, to Perform a difference operation to obtain the termination threshold of the new window loop. That is, at least three window parameter estimates should be performed, when two consecutive... When the difference is less than the preset threshold, it indicates that the degree of change in the difference between two adjacent parameter estimation results meets the expectations, which means that the parameter estimation results of the tropospheric delay model have become stable at this time.

8. The atmospheric delay correction method for time-series InSAR monitoring data according to claim 1, characterized in that, In step (8), after subtracting the atmospheric delay phase from the initial wrapped phase before atmospheric delay correction, 3D phase unwrapping is performed again, and the change in the percentage of error-free points between the two unwrapping operations relative to the total number of interferogram measurement points is calculated as ΔU. Among them, U i+1 and U i N represents the percentage of error-free unwrapped points in all interferograms after the (i+1)th and ithth estimations of the atmospheric delay phase and the removal of the interference phase, respectively; mp This indicates the number of measurement points; M is the number of interferograms; error-free points refer to points where the absolute value of the unwrapped closed-loop residual is less than 1 rad.

9. A computer system, characterized in that, include: One or more processors; A memory that stores operable instructions; The instructions, when executed by the one or more processors, implement the atmospheric delay correction method for time-series InSAR monitoring data according to any one of claims 1-8.

10. A computer-readable medium for storing computer programs, characterized in that, The computer program includes instructions executable by one or more computers, which, when executed by the one or more computers, complete the process of the atmospheric delay correction method for time-series InSAR monitoring data according to any one of claims 1-8.

Images