A time series InSAR atmospheric delay correction method based on quadtree and joint model

By using quadtree segmentation and joint model estimation to estimate atmospheric delay, deformation, and terrain errors, the problem of atmospheric delay error in InSAR was solved, and high-precision deformation monitoring was achieved.

CN116148783BActive Publication Date: 2026-04-24TONGJI UNIV
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
TONGJI UNIV
Filing Date
2022-12-31
Publication Date
2026-04-24

AI Technical Summary

Technical Problem

In existing InSAR technology, atmospheric delay errors can mask the true deformation in areas with large topographic relief, leading to changes in phase gradient and affecting monitoring accuracy and reliability. Furthermore, existing methods struggle to effectively separate atmospheric delay and topographic error signals.

Method used

A quadtree is used to segment the interferogram, and a joint model is used to jointly model and estimate atmospheric delay, deformation and topographic error. Atmospheric delay is corrected by estimating the polynomial parameters in the quadtree sub-window and unwrapping the phase of the sparse linear equation.

Benefits of technology

It improves the accuracy of InSAR large-scale deformation monitoring, is applicable to various deformation scenarios and complex meteorological conditions, and enhances the accuracy and consistency of atmospheric delay correction.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN116148783B_ABST
    Figure CN116148783B_ABST
Patent Text Reader

Abstract

The application relates to a time-series InSAR atmospheric delay correction method based on a quadtree and a joint model, which comprises the following steps: S1, using a phase standard deviation to divide an interferogram into quadtree windows; S2, in each quadtree sub-window after division, jointly modeling and estimating atmospheric delay, deformation and terrain error; and S3, in each quadtree sub-window, based on estimated atmospheric delay polynomial parameters, correcting atmospheric delay phases of each interferogram. Compared with the prior art, the application improves the model coupling degree between atmospheric delay signals and terrain by window division, constructs a joint model according to the time and space characteristics of different signals to synchronously estimate various parameters, thereby establishing an adaptive atmospheric delay correction method and improving the InSAR large-range deformation monitoring precision.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the application of synthetic aperture radar interferometry technology in time-series data processing, and in particular to a time-series InSAR atmospheric delay correction method based on quadtrees and joint models. Background Technology

[0002] Synthetic Aperture Radar Interferometry (InSAR), a novel space-based geodesy technique, features all-weather, wide-area coverage, and high spatial resolution. It retrieves surface deformation information within a time interval by comparing the phase difference between radar images acquired at different times covering the same area. In recent years, this technique has been widely used to monitor various surface deformations caused by natural factors (such as earthquakes, volcanoes, and landslides) or human factors (such as groundwater extraction, mining, and underground tunnel excavation). However, the phase difference obtained by differentially analyzing two radar images in this technique not only contains deformation information but is also affected by topographic errors, atmospheric delay, orbital errors, and decoherence noise. Among these, atmospheric delay, as a major source of error, produces significant phase gradient changes in areas with large topographic relief, thus masking the true deformation and severely limiting the widespread application and reliability of InSAR deformation monitoring.

[0003] Currently, methods for separating and correcting InSAR atmospheric delay errors are mainly divided into two categories: one category uses external data (such as GPS, satellite spectrometers, and numerical atmospheric models) to simulate the atmospheric delay at the time of synthetic aperture radar (SAR) data acquisition. However, this type of method depends on the accuracy of the externally acquired data and is limited by the accuracy of interpolation of discrete data in the time and spatial domains. The other category uses InSAR phase observations itself and sets up a spatiotemporal three-dimensional filter based on prior assumptions to separate atmospheric delay signals. However, this type of method is limited by the universality of the prior assumptions and the human experience setting of the filter window size.

[0004] InSAR atmospheric delay phase is primarily caused by spatiotemporal disturbances in the tropospheric refractive index. Physically, atmospheric delay is influenced by tropospheric water vapor and its spatial distribution is correlated with topography. However, due to differences in weather and climate conditions across regions, the correlation between atmospheric delay phase and topography exhibits non-uniformity across a large spatial range, resulting in non-unique linear regression coefficients at different spatial locations, with some regions even showing diametrically opposed regression coefficients. Furthermore, since InSAR interferometric phase itself contains multiple signal sources such as deformation, atmospheric delay, and topographic errors, relying solely on topographic correlation is insufficient to effectively separate the phase contributions of different signal sources. Therefore, overcoming the spatially significant differences in the regression coefficients of atmospheric delay and topography, and offsetting the bias contributions of other signal sources to the multivariate regression, are key challenges for achieving high-precision InSAR atmospheric delay correction. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the existing technology and provide a temporal InSAR atmospheric delay correction method based on quadtree and joint model.

[0006] The objective of this invention can be achieved through the following technical solutions:

[0007] A temporal InSAR atmospheric delay correction method based on quadtree and joint model includes the following steps:

[0008] S1. Use the phase standard deviation to divide the interferogram into a quadtree window;

[0009] S2. Within each quadtree sub-window after segmentation, atmospheric delay, deformation, and terrain error are jointly modeled and estimated;

[0010] S3. Within each quadtree sub-window, the atmospheric delay phase of each interferogram is corrected based on the estimated atmospheric delay polynomial parameters.

[0011] Furthermore, the specific division of step S1 includes:

[0012] S101. The interferogram is divided into four quadtree sub-windows with overlapping regions using the phase standard deviation.

[0013] S102. Calculate the average phase standard deviation of the M interferograms within each quadtree sub-window.

[0014] S103, If the average phase standard deviation within the quadtree sub-window Greater than the set threshold S threshold Then repeat steps S101 and S102 to further divide the sub-window into four sub-windows with overlapping regions. If the average phase standard deviation within the window... Less than the set threshold S threshold If so, then stop splitting the window.

[0015] Furthermore, in step S2, the atmospheric delay signal is represented in the interferogram as the differential phase of the atmospheric delay at the two SAR image acquisition times, and is modeled in each SAR image using a spatial polynomial:

[0016]

[0017] Where i is the SAR image sequence number. The atmospheric delay phase at the time of image acquisition is given, where X and y represent the range and azimuth coordinates in the SAR image coordinate system, respectively, H is the elevation, and a i b i c i di as well as Let be the coefficients of the polynomial to be determined;

[0018] For M interferograms, the contribution of all atmospheric delay phases within a single quadtree subwindow is expressed as:

[0019]

[0020] In the formula, This represents the Kronecker tensor product. Since InSAR phase is a relative quantity of observation, assuming that the atmospheric delay polynomial parameter corresponding to one of the reference SAR images is known to be 0, P trop Let F be the vector of atmospheric delay polynomial parameters to be determined for the remaining N SAR images, and let F be the transformation matrix between the SAR image and the interferogram, i.e., relative to the reference image.

[0021]

[0022] In this matrix, each row represents an interferogram, and each column represents a SAR image; where 1 and -1 represent the main image and sub-image in the interferogram, respectively. Since the atmospheric delay polynomial parameter of the reference image is known, the column corresponding to the reference image in the F matrix is ​​also deleted.

[0023] Furthermore, in step S2, the deformation signal is modeled based on the correlation between the deformation signal and the time dimension, and the deformation dynamics model is represented by a cubic polynomial.

[0024] Furthermore, in step S2, the terrain error is modeled based on the linear relationship between the terrain error contribution phase and the vertical baseline of the interferogram. The deformation of a pixel in the quadtree sub-window and the terrain error phase contribution are represented as follows:

[0025]

[0026] In the formula, p is the pixel sequence number, and i and j represent the main and sub-images of the interferogram. and These correspond to the three unknown parameters to be determined in the dynamic model: velocity, acceleration, and rate of change of acceleration, respectively. λ is the radar microwave wavelength, r is the distance from the SAR satellite to the ground, and θ is the satellite incident angle. For the vertical baseline of the interferogram, Δh p For unknown terrain errors; considering M interferograms, the formula is extended to:

[0027] φ defo+top,p =[G·D defo D topo ]·P defo+topo,p

[0028] in,

[0029]

[0030]

[0031]

[0032] In the formula, P defo+topo,p represents the unknown deformation and topographic error parameters, G is the conversion matrix between the SAR image interval and the interferogram. If there are S pixels in a quadtree sub-window, the deformation and topographic error phase contributions of all pixels in this quadtree sub-window can be expressed as:

[0033]

[0034] In the formula, I S×S is the identity matrix.

[0035] Further, the conversion matrix G between the SAR image interval and the interferogram is specifically expressed as: for the matrix element G(k, l), if the k-th interferogram is composed of SAR images i and j (i < j), for the elements greater than or equal to i and less than or equal to j, then G(k, l) = t l+1 -t l , otherwise G(k, l) = 0, where t l represents the acquisition date of the l-th SAR image.

[0036] Further, when considering the atmospheric delay, deformation and topographic signals in the quadtree sub-window at the same time, for M interferograms and the unwrapped phases of S pixels, the phase observation model adopts a large-scale sparse linear equation, and the unknown parameters are obtained by Golub-Kahan bidiagonalization least squares, which is specifically expressed as:

[0037]

[0038] In the formula, Φ is the known phase observation, [B defo+topo , B [[ID=四十二]] trop is the known design matrix, which is obtained by combining the design matrix of the atmospheric delay polynomial model and the deformation and topographic error models, is the unknown vector to be solved, including the atmospheric delay polynomial parameters, the deformation dynamics model parameters and the topographic error parameters.

[0039] Further, in step S3, within each quadtree sub-window, the atmospheric delay phase of each interferogram is corrected based on the estimated atmospheric delay polynomial parameters, that is:

[0040]

[0041] In the formula, The original phase of the interferogram. This is the phase of the interferogram after atmospheric delay correction.

[0042] Furthermore, after the atmospheric delay phase within each quadtree sub-window is corrected, the phases of multiple quadtree sub-windows are merged and stitched together to form the overall interferogram phase. A Delaunay triangulation is established within each quadtree sub-window, and the phase difference on the arc segments of the triangulation is calculated. Since there is pixel overlap between adjacent quadtree sub-windows, the triangulations constructed by adjacent windows will also overlap, meaning there will be common arc segments within adjacent windows. Considering that independent parameter estimation within each quadtree sub-window will lead to inconsistent phase difference values ​​for common arc segments, the average value of the phase difference values ​​of common arc segments in different windows can be used to replace the common arc segments. Then, the atmospheric delay-corrected interferogram phase is obtained by least squares.

[0043] Furthermore, the phase difference between the atmospheric delay-corrected phase of the entire interferogram and the phase difference on the arc segment is specifically expressed as follows:

[0044]

[0045] In the formula, For the observation of phase difference values ​​on the arc segment, Let C be the phase of the interferogram after atmospheric delay correction. Let C be the transformation matrix between pixels and arcs. The number of rows in matrix C is equal to the total number of arcs, and the number of columns is equal to the total number of pixels in the interferogram minus one. The start and end points of each arc are represented by 1 and -1, respectively, and the remaining matrix elements are 0. It should be noted that, in order to avoid rank deficiency in matrix C, the column corresponding to the reference point has been deleted.

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

[0047] 1. This invention improves the model coupling between atmospheric delay signals and terrain by using quadtree segmentation interferograms, and constructs a joint model to synchronously estimate various parameters by utilizing the spatiotemporal characteristics of atmospheric delay, deformation and terrain error components, thereby establishing an adaptive atmospheric delay correction method to improve the accuracy of InSAR large-scale deformation monitoring.

[0048] 2. This invention uses the phase standard deviation of interferograms as an indicator and employs a quadtree to segment the interferogram sheet, making full use of the local characteristics of tropospheric atmospheric delay and terrain correlation, thereby improving the model coupling between atmospheric delay signals and terrain.

[0049] 3. This invention utilizes the spatiotemporal statistical characteristics of time-series phase to establish a joint model for synchronous estimation of atmospheric delay, deformation, and topographic error parameters, thereby improving the estimation accuracy of atmospheric delay polynomial parameters. This allows the atmospheric delay correction method to be applied to a wide range of deformation scenarios and areas with complex meteorological conditions, providing high-precision data for deformation monitoring. Attached Figure Description

[0050] Figure 1 This is a flowchart of the temporal InSAR atmospheric delay correction of the present invention;

[0051] Figure 2 This is an example diagram of the quadtree segmentation of the interferogram image according to the present invention;

[0052] Figure 3 This is a comparison of the interferogram before and after atmospheric delay correction according to the present invention, wherein the first row is the original interferogram and the second row is the interferogram after atmospheric delay correction;

[0053] Figure 4 This is a statistical graph of the phase standard deviation of the interferogram before and after atmospheric delay correction in this invention. Detailed Implementation

[0054] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. These embodiments are based on the technical solution of the present invention and provide detailed implementation methods and specific operating procedures. However, the scope of protection of the present invention is not limited to the following embodiments.

[0055] Example

[0056] like Figure 1 As shown, a time-series InSAR atmospheric delay correction method based on quadtree and joint model includes the following steps:

[0057] S1. Use the phase standard deviation to divide the interferogram into a quadtree window;

[0058] S2. Within each quadtree sub-window after segmentation, atmospheric delay, deformation, and terrain error are jointly modeled and estimated;

[0059] S3. Within each quadtree sub-window, the atmospheric delay phase of each interferogram is corrected based on the estimated atmospheric delay polynomial parameters.

[0060] The specific segmentation in step S1 includes:

[0061] S101. The interferogram is divided into four quadtree sub-windows with overlapping regions using the phase standard deviation.

[0062] S102. Calculate the average phase standard deviation of the M interferograms within each quadtree sub-window.

[0063] S103, If the average phase standard deviation within the quadtree sub-window Greater than the set threshold S threshold Then repeat steps S101 and S102 to further divide the sub-window into four sub-windows with overlapping regions. If the average phase standard deviation within the window... Less than the set threshold S threshold If so, then stop splitting the window.

[0064] In step S2, the atmospheric delay signal can be represented in the interferogram as the differential phase of the atmospheric delay at the two SAR image acquisition times, and modeled in each SAR image using a spatial polynomial:

[0065]

[0066] Where i is the SAR image sequence number. The atmospheric delay phase at the time of image acquisition is given, where X and Y represent the range and azimuth coordinates in the SAR image coordinate system, respectively, H is the elevation, and a i b i c i d i as well as Let be the coefficients of the polynomial to be determined;

[0067] For M interferograms, the contribution of all atmospheric delay phases within a single quadtree subwindow is expressed as:

[0068]

[0069] In the formula, This represents the Kronecker tensor product. Since InSAR phase is a relative quantity of observation, assuming that the atmospheric delay polynomial parameter corresponding to one of the reference SAR images is known to be 0, P trop Let F be the vector of atmospheric delay polynomial parameters to be determined for the remaining N SAR images, and let F be the transformation matrix between the SAR image and the interferogram, i.e., relative to the reference image.

[0070]

[0071] In this matrix, each row represents an interferogram, and each column represents a SAR image; where 1 and -1 represent the main image and sub-image in the interferogram, respectively. Since the atmospheric delay polynomial parameter of the reference image is known, the column corresponding to the reference image in the F matrix is ​​also deleted.

[0072] In step S2, the deformation signal is modeled based on its correlation with the time dimension, and the deformation dynamics model is represented by a cubic polynomial.

[0073] In step S2, for the terrain error, a model is established based on the linear relationship between the terrain error contribution phase and the perpendicular baseline of the interferogram. The deformation and terrain error phase contribution of a pixel in the quadtree sub-window are expressed as:

[0074]

[0075] In the formula, p is the pixel serial number, i and j represent the master and slave images of the interferogram, and correspond to the three unknown parameters to be determined, namely the rate, acceleration, and acceleration change rate in the dynamic model, λ is the radar microwave wavelength, r is the distance from the SAR satellite to the ground, θ is the satellite incident angle, is the perpendicular baseline of the interferogram, and Δh p is the unknown terrain error; considering M interferograms, the formula is extended to:

[0076] φ defo+topo,p =[G·D defo ,D topo ·P defo+topo,p

[0077] Among them,

[0078]

[0079]

[0080]

[0081] In the formula, P defo+topo,p represents the unknown deformation and terrain error parameters, G is the conversion matrix between the SAR image interval and the interferogram. If there are S pixels in a quadtree sub-window, the deformation and terrain error phase contributions of all pixels in this quadtree sub-window can be expressed as:

[0082]

[0083] In the formula, I S×S is the identity matrix.

[0084] The conversion matrix G between the SAR image interval and the interferogram is specifically expressed as: for the matrix element G(k, l), if the kth interferogram is composed of SAR images i and j (i < j), for elements greater than or equal to i and less than or equal to j, then G(k, l) = t l+1 -t l , otherwise G(k, l) = 0, where t l represents the acquisition date of the lth SAR image.

[0085] When atmospheric delay, deformation, and topographic signals are considered simultaneously within a quadtree sub-window, for M interferograms and S unwrapped phases, the phase observation model employs a large-scale sparse linear equation. The unknown parameters are obtained through Golub-Kahan bidiagonalized least squares, specifically expressed as follows:

[0086]

[0087] In the formula, Φ represents the known phase observation, [B defo+top B trop Given the design matrix, it is obtained by combining the design matrix of the atmospheric delay polynomial model and the deformation and terrain error models. The unknown vectors to be determined include atmospheric delay polynomial parameters, deformation dynamics model parameters, and terrain error parameters.

[0088] In step S3, within each quadtree sub-window, the atmospheric delay phase of each interferogram is corrected based on the estimated atmospheric delay polynomial parameters, i.e.:

[0089]

[0090] In the formula, The original phase of the interferogram. This is the phase of the interferogram after atmospheric delay correction.

[0091] After atmospheric delay phase correction is completed within each quadtree sub-window, the phases of multiple quadtree sub-windows are merged and stitched together to form the overall interferogram phase. A Delaunay triangulation is built within each quadtree sub-window, and the phase difference on the arc segments of the triangulation is calculated. Since there is pixel overlap between adjacent quadtree sub-windows, the triangulations constructed by adjacent windows will also overlap, meaning there will be common arc segments within adjacent windows. Considering that independent parameter estimation within each quadtree window will lead to inconsistent phase difference values ​​for common arc segments, the average phase difference value of common arc segments in different windows can be used as a substitute. Then, the phase of the interferogram after atmospheric delay correction is obtained by least squares.

[0092] The phase difference between the atmospheric delay correction phase and the phase on the arc segment of the entire interferogram is specifically expressed as follows:

[0093]

[0094] In the formula, For the observation of phase difference values ​​on the arc segment, Let C be the phase of the interferogram after atmospheric delay correction. Let C be the transformation matrix between pixels and arcs. The number of rows in matrix C is equal to the total number of arcs, and the number of columns is equal to the total number of pixels in the interferogram minus one. The start and end points of each arc are represented by 1 and -1, respectively, and the remaining matrix elements are 0. It should be noted that, in order to avoid rank deficiency in matrix C, the column corresponding to the reference point has been deleted.

[0095] The experiment used 54 Sentinel-1A ascending orbit image data acquired from the Ete Are volcano region in Ethiopia, Africa, spanning from February 4, 2017 to September 4, 2019. After preprocessing, 53 differential interferometric unwrapped phase maps were generated. These interferograms were then segmented using a quadtree, and the results are shown below. Figure 2 As shown, the interferogram phase atmospheric delay correction is compared before and after, for example... Figure 3 As shown, the statistical values ​​of the phase standard deviation of the interferograms before and after atmospheric delay correction are as follows: Figure 4 As shown. Experimental results show that the present invention can effectively correct the atmospheric delay phase in time-series InSAR interferograms, reducing the average standard deviation of the interferogram phase from 4.6 radians before correction to 1.3 radians after correction.

[0096] The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and variations based on the concept of the present invention without creative effort. Therefore, all technical solutions that can be obtained by those skilled in the art based on the concept of the present invention through logical analysis, reasoning, or limited experimentation on the basis of existing technology should be within the scope of protection defined by the claims.

Claims

1. A temporal InSAR atmospheric delay correction method based on quadtree and joint model, characterized in that, Includes the following steps: S1. Use the phase standard deviation to divide the interferogram into a quadtree window; The specific segmentation in step S1 includes: S101. The interferogram is divided into four quadtree sub-windows with overlapping regions using the phase standard deviation. S102. Calculate the average phase standard deviation of the M interferograms within each quadtree sub-window. ; S103, If the average phase standard deviation within the quadtree sub-window Greater than the set threshold Then repeat steps S101 and S102 to further divide the sub-window into four sub-windows with overlapping regions. If the average phase standard deviation within the window... Less than the set threshold If so, then stop splitting the window; S2. Within each quadtree sub-window after segmentation, atmospheric delay, deformation, and terrain error are jointly modeled and estimated; S3. Within each quadtree sub-window, the atmospheric delay phase of each interferogram is corrected based on the estimated atmospheric delay polynomial parameters.

2. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 1, characterized in that, In step S2, the atmospheric delay signal is represented in the interferogram as the differential phase of the atmospheric delay at the two SAR image acquisition times, and is modeled using a spatial polynomial in each SAR image: in, This is the SAR image sequence number. To obtain the temporal atmospheric delay phase for this image, and These represent the range and azimuth coordinates in the SAR image coordinate system, respectively. For elevation, , , , as well as Let be the coefficients of the polynomial to be determined; For M interferograms, the contribution of all atmospheric delay phases within a single quadtree subwindow is expressed as: In the formula, Let Kronecker tensor product be used. Since InSAR phase is a relative quantity of observation, it is assumed that the atmospheric delay polynomial parameter of one of the reference SAR images is known to be 0. Let be the vector of atmospheric delay polynomial parameters to be determined for the remaining N SAR images. This is the transformation matrix between SAR images and interferograms, i.e., relative to the reference image. In this matrix, each row represents an interferogram, and each column represents a SAR image; where 1 and -1 represent the main image and sub-image in the interferogram, respectively. Since the atmospheric delay polynomial parameter of the reference image is known, The column corresponding to the reference image in the matrix was also deleted.

3. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 1, characterized in that, In step S2, the deformation signal is modeled based on the correlation between the deformation signal and the time dimension, and the deformation dynamics model is represented by a cubic polynomial.

4. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 2, characterized in that, In step S2, the terrain error is modeled based on the linear relationship between the terrain error contribution phase and the vertical baseline of the interferogram. The deformation of a pixel in a quadtree sub-window and the terrain error phase contribution are represented as follows: In the formula, For the pixel sequence number, and This represents the main and secondary images in the interferogram. , and These correspond to the three unknown parameters to be determined in the dynamic model: velocity, acceleration, and rate of change of acceleration. For radar microwave wavelength, This refers to the distance from the SAR satellite to the ground. The angle of incidence of the satellite. The vertical baseline of the interferogram. For unknown terrain errors; consider Aspect ratio interferogram, the formula is extended to: in, In the formula, This represents unknown deformation and terrain error parameters. This is the transformation matrix between SAR image intervals and interferograms. If a quadtree sub-window contains... If there are 100 pixels, then the deformation and terrain error phase contribution of all pixels within the quadtree sub-window can be expressed as: In the formula, It is an identity matrix.

5. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 4, characterized in that, The transformation matrix between the SAR image interval and the interferogram Specifically, for matrix elements If the first Interferograms from SAR images and composition, For greater than or equal to and less than or equal to The element, then ,otherwise ,in Indicates the first The date of acquisition of the SAR image.

6. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 4, characterized in that, When atmospheric delay, deformation, and terrain signals are considered simultaneously within the quadtree sub-window, for M interferograms and The phase observation model for each pixel unwrapped phase employs a large sparse linear equation, and the unknown parameters are obtained through Golub-Kahan bidiagonalized least squares, specifically expressed as follows: In the formula, For observations with known phase, Given the design matrix, it is obtained by combining the design matrix of the atmospheric delay polynomial model and the deformation and terrain error models. The unknown vectors to be determined include atmospheric delay polynomial parameters, deformation dynamics model parameters, and terrain error parameters.

7. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 4, characterized in that, In step S3, within each quadtree sub-window, the atmospheric delay phase of each interferogram is corrected based on the estimated atmospheric delay polynomial parameters, i.e.: In the formula, The original phase of the interferogram. The phase of the interferogram after atmospheric delay correction.

8. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 7, characterized in that, After the atmospheric delay phase within each quadtree sub-window is corrected, the phases of multiple quadtree sub-windows are merged and stitched together to form the overall interferogram phase. A Delaunay triangulation is established within each quadtree sub-window, and the phase difference on the arc segments of the triangulation is calculated. Since there is pixel overlap between adjacent quadtree sub-windows, the triangulations constructed by adjacent windows will also overlap, meaning there will be common arc segments within adjacent windows. Considering that independent parameter estimation within each quadtree window will lead to inconsistent phase difference values ​​for common arc segments, the average phase difference value of common arc segments in different windows can be used as a substitute. Then, the atmospheric delay-corrected interferogram phase is obtained by least squares.

9. The temporal InSAR atmospheric delay correction method based on quadtree and joint model according to claim 8, characterized in that, The phase difference between the atmospheric delay-corrected phase and the phase on the arc segment of the entire interferogram is specifically expressed as follows: In the formula, For the observation of phase difference values ​​on the arc segment, To determine the phase of the interferogram after atmospheric delay correction, Let be the transformation matrix between pixels and arc segments, where The number of rows in the matrix equals the total number of arc segments, and the number of columns equals the total number of pixels in the interferogram minus one. The start and end points of each arc segment are represented by 1 and -1, respectively, and the remaining matrix elements are 0. It is important to note that, to avoid... Rank deficiency of a matrix The column corresponding to the reference point in the matrix has been deleted.