A method for error estimation of InSAR time series DEM

Abstract

The present invention discloses an InSAR time series DEM error estimation method, which belongs to the field of terrain monitoring technology and can solve the deformation time series and deformation rate errors caused by DEM errors in different time periods in traditional SBAS-InSAR technology, as well as the problem that the landfill height in different time periods cannot be obtained. The method comprises: S1, obtaining multiple single-view complex images of the area to be monitored, and performing differential interference processing and phase optimization processing on the multiple single-view complex images to obtain multiple optimized interference graphs; S2, performing phase unwrapping on the multiple interference graphs to obtain multiple unwrapped graphs; S3, using an inversion function to invert the multiple unwrapped graphs into multiple continuous phase graphs, and obtaining n grouping time points of the deformation phase time series of each pixel according to the multiple continuous phase graphs; S4, using the n grouping time points to divide the deformation phase time series of the corresponding pixels into n+1 groups, and obtaining the DEM error corresponding to each group of pixels. The present invention is used for terrain monitoring.

Description

Technical Field

[0001] The invention relates to an InSAR time series DEM error estimation method, belonging to the technical field of terrain monitoring. Background Art

[0002] When solving the deformation time series of the monitoring area, the traditional SBAS-InSAR (Small BAseline Subset-InSAR) technology assumes that there will be no obvious terrain fluctuations on the surface (for example, digging and landfilling), that is, the DEM (Digital Elevation Model) error of the monitoring area is fixed. However, for large-scale projects such as land reclamation and mountain digging, the ground landfill is not completed in one go. After the initial landfill is completed, the consolidation of the loose soil on the seabed and the initial landfill material will cause ground subsidence and cause a DEM error Δh1. After a period of time after the initial landfill is completed, another landfill will be carried out. At this time, the consolidation of the secondary landfill material increases the pressure on the previous landfill material and the loose soil on the seabed, which will cause the acceleration of ground subsidence and cause a new DEM error Δh2. Therefore, after multiple landfills, multiple different DEM errors Δh1, Δh2, ..., Δh n .

[0003] Traditional SBAS-InSAR technology can only estimate a DEM error ΔH, but for land reclamation projects that involve multiple landfills, there are different DEM errors in different time periods, and the landfill time is also different for points at different ground locations. Therefore, the use of traditional SBAS-InSAR technology will affect the calculation accuracy of the deformation time series and deformation rate of the monitored area. At the same time, when solving such problems, traditional technology cannot solve the landfill height and landfill process in different time periods, nor can it quantitatively explore the spatiotemporal landfill process of each landfill project. Summary of the invention

[0004] The present invention provides an InSAR time series DEM error estimation method, which can solve the problems of deformation time series and deformation rate errors caused by DEM errors in different time periods in traditional SBAS-InSAR technology, as well as the inability to obtain landfill heights in different time periods.

[0005] The present invention provides an InSAR time series DEM error estimation method, the method comprising:

[0006] S1, acquiring a plurality of single-view complex images of the area to be monitored, and performing differential interference processing and phase optimization processing on the plurality of single-view complex images to obtain a plurality of optimized interference patterns;

[0007] S2, performing phase unwrapping on the multiple interference patterns to obtain multiple unwrapping patterns;

[0008] S3, inverting the plurality of unwrapped images into a plurality of continuous phase images using an inversion function, and obtaining n grouped time points of the deformation phase time series of each pixel according to the plurality of continuous phase images;

[0009] S4. Divide the deformation phase time series of the corresponding pixels into n+1 groups using the n grouping time points, and obtain the DEM error corresponding to each group of pixels.

[0010] Optionally, the step S1 performs differential interference processing and phase optimization processing on the plurality of single-view complex images to obtain a plurality of optimized interference patterns, specifically including:

[0011] S11, preprocessing the plurality of single-view complex images, and performing differential interference processing on the preprocessed images in pairs to obtain a plurality of interference patterns;

[0012] S12, using SRTM-DEM to remove the external terrain phase in the interference pattern to obtain a processed interference pattern;

[0013] S13, extracting homogeneous points in the processed interference pattern, and using the homogeneous points to perform phase noise reduction on the interference pattern to obtain a plurality of optimized interference patterns.

[0014] Optionally, the S2 specifically includes:

[0015] S21, performing MCF phase unwrapping on the optimized multiple interference patterns to obtain multiple unwrapped patterns;

[0016] S22, performing unwrapping error identification and phase correction on each pixel in the multiple unwrapping images to obtain multiple unwrapping images after phase correction.

[0017] Optionally, the S22 specifically includes:

[0018] S221, normalizing corresponding pixels in a plurality of unwrapped images to obtain deformation phase rates of the pixels;

[0019] S222, sorting the deformation phase rates of the pixels according to the middle time points of the corresponding interference patterns, and obtaining the median rate of the deformation phase rates of the sorted pixels in a time window of a preset time length before and after the current moment;

[0020] S223, using the median rate to perform phase correction on the deformation phase rate of the pixel at the current moment, to obtain a plurality of phase-corrected unwrapping images.

[0021] Optionally, the S3 specifically includes:

[0022] S31, selecting an unwrapping graph that meets preset requirements from multiple unwrapping graphs after phase correction;

[0023] S32, using an inversion function to invert the unwrapped image that meets the preset requirements into a plurality of continuous phase images;

[0024] S33, performing denoising on the multiple continuous phase images using a one-dimensional wavelet transform method to obtain a denoised deformation phase time series of each pixel;

[0025] S34. Obtain n grouping time points according to the deformation phase time series of each pixel.

[0026] Optionally, the S34 is specifically:

[0027] The first-order derivative of the deformation phase time series of each pixel is taken to obtain n peak values ​​in the deformation phase time series of each pixel, and the time points corresponding to the peak values ​​are the grouping time points.

[0028] Optionally, the S34 is specifically:

[0029] S341, determining a sliding window from the deformation phase time series of each pixel, and fitting the deformation phase within the sliding window to obtain a fitting function;

[0030] S342, predicting the deformation phase of a future time period adjacent to the sliding window according to the fitting function to obtain a predicted phase time series;

[0031] S343: Obtain the measured phase time series of the future time period, calculate the standard deviation or variance between the measured phase time series and the predicted phase time series, and use the peak value of the standard deviation or variance as the grouping time point.

[0032] Optionally, the S4 specifically includes:

[0033] S41, using the n grouping time points to divide the deformation phase time series of the corresponding pixels into n+1 groups;

[0034] S42, selecting the unwrapped phases that meet the preset conditions in each group of pixels, constructing a DEM solution equation for the corresponding pixel group according to the unwrapped phases that meet the preset conditions, and obtaining the DEM error corresponding to each group of pixels according to the DEM solution equation.

[0035] Optionally, after S4, the method further includes:

[0036] S5. De-noise the DEM errors corresponding to each group of pixels, and generate a ground elevation landfill time series based on the de-noised DEM errors.

[0037] Optionally, after S4, the method further includes:

[0038] S6, denoising the DEM error corresponding to each group of pixels, and removing the corresponding denoised DEM error from the unwrapped phase of the unwrapped image to obtain an unwrapped image containing only the deformation phase;

[0039] S7. Obtaining the deformation time series and deformation rate of the monitoring area according to the unwrapped image containing only the deformation phase.

[0040] Optionally, the S223 specifically includes:

[0041] Use the first formula to perform phase correction;

[0042] The first formula is:

[0043]

[0044] Among them, UNW t Indicates the corrected phase; fix() indicates rounding to 0; TM t The length of the time base representing the phase at time t; is the median rate; is the deformation phase rate at time t.

[0045] Optionally, the S7 is specifically:

[0046] According to the unwrapping graph containing only the deformation phase, the deformation time series of the monitoring area is obtained by using the least square method or the SVD decomposition method, and the deformation rate of the monitoring area is obtained by using the phase superposition method.

[0047] The beneficial effects that the present invention can produce include:

[0048] The InSAR time series DEM error estimation method provided by the present invention can not only obtain the deformation time series and deformation rate of the monitoring area, but also obtain the elevation time series of the time series ground landfill, and use the obtained elevation deformation time series to correct the DEM error in the deformation time series, thereby improving the accuracy of deformation monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 A flow chart of the InSAR time series DEM error estimation method provided by an embodiment of the present invention;

[0050] Figure 2 A block diagram of an InSAR time series DEM error estimation method provided by an embodiment of the present invention;

[0051] Figure 3A time-space baseline diagram in the experimental results provided by an embodiment of the present invention;

[0052] Figure 4 A schematic diagram of the cumulative DEM error in the experimental results provided by the embodiment of the present invention;

[0053] Figure 5 The DEM error time series in the experimental results provided by the embodiment of the present invention;

[0054] Figure 6 The DEM error deformation time series of the main regions in the experimental results provided by the embodiment of the present invention;

[0055] Figure 7 for Figure 4 Comparison chart of DEM error correction and original time series at point D3;

[0056] Figure 8 for Figure 4 Comparison chart of DEM error correction and original time series at point D4. DETAILED DESCRIPTION

[0057] The present invention is described in detail below in conjunction with embodiments, but the present invention is not limited to these embodiments.

[0058] The embodiment of the present invention provides an InSAR time series DEM error estimation method, such as Figure 1 and Figure 2 As shown, the method includes:

[0059] S1. Acquire multiple single-view complex images of the area to be monitored, and perform differential interference processing and phase optimization processing on the multiple single-view complex images to obtain multiple optimized interference patterns.

[0060] After acquiring multiple single-view complex images of the area to be monitored, S1 specifically includes:

[0061] S11, preprocessing a plurality of single-view complex images, and combining the preprocessed images in pairs to perform differential interference processing to obtain a plurality of interference patterns.

[0062] There are many ways of preprocessing, and those skilled in the art can set them according to actual conditions, which are not limited in the embodiments of the present invention. In practical applications, the Sentinel-1 data (i.e., single-view complex SLC images) are preprocessed, including radiation correction, burst cropping and splicing, orbit refinement using precise orbit data, and alignment to 1 / 1000 pixel and azimuth spectrum de-skew with the assistance of external SRTM (Shuttle Radar Topography Mission).

[0063] After the preprocessing is completed, the preprocessed images are combined in pairs, and differential interferometry processing is performed using a multi-view ratio of 4:1 in the range and azimuth directions to obtain multiple interferograms. At this time, the corresponding surface resolution is approximately 16mx16m.

[0064] S12. Use SRTM-DEM to remove the external terrain phase in the interference pattern to obtain a processed interference pattern.

[0065] S13, extracting homogeneous points in the processed interference pattern, and using the homogeneous points to perform phase noise reduction on the interference pattern to obtain multiple optimized interference patterns.

[0066] For DS (distributed scatterer) targets in the processed interferogram, the Anderson-Darling (AD) test method can be used to extract SHPS (homogeneous points, i.e., targets of the same type). After that, phase optimization can be performed by decomposing the robustly estimated covariance matrix using homogeneous points. Specifically, the normalized complex coherence matrix T is an orthogonal matrix that can be decomposed into:

[0067]

[0068] Where N represents the velocity of SLC and H represents the transpose of the matrix. i is a non-negative eigenvalue, μ i is the corresponding eigenvector; without loss of generality, let λ1>λ2>…>λ N ,λ i The larger the value, the more dominant the scattering phase is. i Corresponding to an independent scattering mechanism, T can be decomposed into the signal phase T signal With noise phase T noise ,Right now

[0069]

[0070] Therefore, after performing eigenvalue decomposition on the robustly estimated covariance matrix, the largest eigenvalue and eigenvector, i.e., the dominant scatterer, are found, thereby achieving the purpose of phase denoising of the differential interferogram.

[0071] S2. Phase unwrapping is performed on the multiple interference patterns to obtain multiple unwrapped patterns.

[0072] Specifically include:

[0073] S21, performing MCF phase unwrapping on the optimized multiple interference patterns to obtain multiple unwrapped patterns.

[0074] Since the deformation of ground subsidence is approximately linear in a very short time, the deformation rate is a constant in a very short time. Therefore, 3D phase unwrapping can be performed on all phase-optimized interference graphs. For example, MCF (Minimum Cost Flow) phase unwrapping can be performed to obtain multiple unwrapped graphs.

[0075] S22, performing unwrapping error identification and phase correction on each pixel in the multiple unwrapping images to obtain multiple unwrapping images after phase correction.

[0076] Specifically include:

[0077] S221, normalize the corresponding pixels in the multiple unwrapped images to obtain the deformation phase rate of the pixels.

[0078] Specifically, for a common pixel (X, Y) in the unwrapped images obtained by all MCF calculations, the unwrapped phase is divided by the corresponding time baseline length to obtain the deformation phase rate of the pixel per day.

[0079] S222, sorting the deformation phase rates of the pixels according to the middle time points of the corresponding interference patterns, and obtaining the median rate of the deformation phase rates of the sorted pixels in a time window of a preset time length before and after the current moment.

[0080] The preset duration is a pre-set time length, which can be set by those skilled in the art according to actual conditions, and is not limited in the embodiments of the present invention. For example, the preset duration can be set to 15 days, 18 days, or 20 days, etc.

[0081] S223, using the median rate to perform phase correction on the deformation phase rate of the pixel at the current moment, to obtain a plurality of phase-corrected unwrapping images.

[0082] In the embodiment of the present invention, the deformation phase rate of the pixel is sorted according to the middle time point sequence of the unwrapping graph. Since the median of a set of data has a good anti-error effect, when judging whether there is an unwrapping error in the unwrapping phase of the point (X, Y) at time t, for example, the median rate of the interference on the unwrapping deformation phase rate within a time window of 18 days before and after time t (i.e., the current time) can be selected. Then the deformation phase rate at time t For comparison, the following first formula is used for phase correction.

[0083] Specifically, the first formula is:

[0084]

[0085] Among them, UNW tIndicates the corrected phase; fix() indicates rounding to 0; TM t The length of the time base representing the phase at time t; is the median rate; is the deformation phase rate at time t.

[0086] S3. Invert the multiple unwrapped images into multiple continuous phase images using an inversion function, and obtain n grouped time points of the deformation phase time series of each pixel according to the multiple continuous phase images.

[0087] Due to the consolidation of the secondary reclamation materials, the pressure on the previous landfill materials and the loose soil on the seabed increases, which will cause the acceleration of ground subsidence. Therefore, this accelerated feature can be used to find the time node of the second landfill, so as to group the unwrapped map and solve the DEM error.

[0088] Specifically include:

[0089] S31. Selecting an unwrapping graph that meets preset requirements from the multiple unwrapping graphs after phase correction.

[0090] After 3D phase unwrapping, M interference pairs (ie, unwrapping graphs) with high quality (eg, low noise) can be manually selected. Afterwards, a second-order polynomial can be used to remove the influence of the trend surface for all selected unwrapping graphs.

[0091] S32, using an inversion function to invert the unwrapping image that meets preset requirements into a plurality of continuous phase images.

[0092] The selected M unwrapped images can be inverted into N continuous phase images. The pixel-by-pixel inversion function is defined as follows:

[0093] δφ=Bφ;

[0094] φ=(B T B) -1 B T δφ;

[0095] Among them, δφ=[UNW1,UNW2,…,UNW M ] T It is represented as a known vector of the unfolded phase values ​​of the M unwrapped graphs. φ=[φ 1 ,φ 2 ,…,φ N ] T The unknown vector represents the phase values ​​of N consecutive phase maps in the time series, and B is the design matrix connecting the disentangled map and the consecutive phase maps.

[0096] S33, using a one-dimensional wavelet transform method to perform denoising on multiple continuous phase images to obtain a deformation phase time series of each pixel after denoising.

[0097] For a certain pixel deformation phase time series contains observation noise. Therefore, it can be expressed as:

[0098] φ (i,j) (t) = S (i,j) (t)+E (i,j) (t), t=1,2,…,N;

[0099] Among them, S (i,j) (t) represents the real signal; E (i,j) (t) represents white noise. Perform wavelet transform on both sides of the above equation:

[0100] WT φ (a,τ)=WT s (a,τ)+WT E (a,τ);

[0101] Among them, WT φ (a,τ) represents the wavelet transform function, a represents the scale, and τ represents the shift.

[0102]

[0103] Where f(t) represents a square integrable function, * represents a complex conjugate, represents the wavelet function.

[0104] After orthogonal wavelet transform, the signal can be removed to the greatest extent. The correlation of E concentrates most of the energy on a few wavelet coefficients with relatively large amplitudes. (i,j) (t) will be distributed on all time axes at all scales after wavelet transform, and the amplitude is not very large. Using this principle, the wavelet coefficients of noise are reduced to the maximum extent at each scale of wavelet transform, and then the processed wavelet coefficients are used to reconstruct the signal, thereby achieving the purpose of suppressing noise.

[0105] S34. Obtain n grouping time points according to the deformation phase time series of each pixel.

[0106] Specifically, the first-order derivative of the deformation phase time series of each pixel may be taken to obtain n peak values ​​in the deformation phase time series of each pixel, and the time point corresponding to the peak value is the grouping time point.

[0107] In the deformation phase time series φ (i,j) After performing one-dimensional wavelet transform to remove noise, we calculate S (i,j) Take the derivative and get its first-order derivative as follows:

[0108] V (i,j)=DIFF(S (i,j) );

[0109] V (i,j) Represents the first derivative of the deformation phase time series, that is, the velocity. DIFF() stands for the first-order inverse.

[0110] When the landfill is completed again, the deformation of the ground settlement will be significantly accelerated, which is reflected in V (i,j) The peak value appears above. Find V (i,j) The peak value in it means that the grouping time point is found, so the unwrapping map can be grouped, the DEM error of each time period can be solved separately, and the DEM error component can be subtracted from the unwrapping map. (i,j) The purpose of wavelet transform is to remove tiny peak values ​​that may lead to misjudgment as breakpoints.

[0111] In another embodiment of the present invention, S34 is specifically:

[0112] S341, determining a sliding window from the deformation phase time series of each pixel, and fitting the deformation phase within the sliding window to obtain a fitting function;

[0113] S342, predicting the deformation phase of a future time period adjacent to the sliding window according to the fitting function to obtain a predicted phase time series;

[0114] S343: Obtain the measured phase time series of the future time period, calculate the standard deviation or variance between the measured phase time series and the predicted phase time series, and use the peak value of the standard deviation or variance as the grouping time point.

[0115] For the deformation phase time series of each pixel, linear fitting is performed by opening a window in time and predicting the deformation phase in the future time period. The predicted value is compared with the measured value to obtain the standard deviation or variance of the predicted time period (i.e., the future time period), and the peak value of the standard deviation or variance is taken as the grouping time point. For example, linear fitting is performed on 1 to 4 phase time series in the deformation phase time series of a certain pixel, and then the fitting function obtained by fitting is used to predict the phase time series of 5-8 time points, and the predicted value is compared with the true value to obtain the standard deviation or variance, and so on, to obtain the standard deviation sequence or variance sequence corresponding to the deformation phase time series of the pixel, and the peak value of the standard deviation sequence or variance sequence is taken as the grouping time point.

[0116] S4. Divide the deformation phase time series of the corresponding pixels into n+1 groups using n grouping time points, and obtain the DEM error corresponding to each group of pixels.

[0117] Specifically include:

[0118] S41, using n grouping time points to divide the deformation phase time series of corresponding pixels into n+1 groups.

[0119] In step S3, the time point at which a pixel (i, j) group is found in the disentangled graph After that, the unwrapped phase can be automatically divided into n+1 groups according to this time node.

[0120] S42, selecting the unwrapped phases that meet the preset conditions in each group of pixels, constructing the DEM solution equations for the corresponding pixel groups according to the unwrapped phases that meet the preset conditions, and obtaining the DEM errors corresponding to each group of pixels according to the DEM solution equations.

[0121] Before constructing the DEM error solution equation group, it is necessary to select the untangling pairs that meet the preset conditions according to certain coherence, time baseline and vertical baseline thresholds. For example, the preset conditions can be set as: coherence>0.6, time baseline≤24 days, vertical baseline>70 meters.

[0122] After selecting the untangle pair, the DEM solution equation can be constructed according to the following formula.

[0123] V w =A w X w -L w , w=1,2,…,n+1;

[0124]

[0125] in,

[0126] V w =[v1,v2,…,v s ] T ;

[0127]

[0128] L w =[UNW1,UNW2,…,UNW s, ] T ;

[0129] X w =[Δh w ,Δv w ] T ;

[0130] Where Δv w , Δh w are the deformation rate and DEM error of the wth group respectively; s is the number of unwrapped pairs that meet the conditions of the wth group; UNW i(i=1,…,s) represents the unwrapped phase; λ, R and θ are the radar wavelength, the distance from the sensor to the target and the incident angle, respectively. and T i (i=1,…,s) are the vertical baseline and time baseline respectively.

[0131] After S4, the method further includes:

[0132] S5. De-noise the DEM errors corresponding to each group of pixels, and generate a ground elevation landfill time series based on the de-noised DEM errors.

[0133] In practical applications, the obtained Δh i , (i=1,2,…,s) median filtering is performed to generate the ground elevation landfill time series.

[0134] Furthermore, after S4, the method further includes:

[0135] S6. De-noise the DEM error corresponding to each group of pixels, and remove the corresponding de-noised DEM error from the unwrapped phase of the unwrapped map to obtain an unwrapped map containing only the deformation phase.

[0136] After calculating the time series DEM error, the following formula can be used to subtract the component related to the DEM error from the unwrapped phase.

[0137]

[0138] in,

[0139]

[0140]

[0141]

[0142] In the above formula It means that the w-th group only contains part of the deformation phase, and g represents the number of all unwrapped phases in the w-th group. It represents the unwrapped phase of the components including DEM error and deformation. w is the corresponding design matrix.

[0143] S7. Obtain the deformation time series and deformation rate of the monitoring area according to the unwrapped map containing only the deformation phase.

[0144] Specifically: Based on the unwrapping graph containing only the deformation phase, the least square method or SVD decomposition (Singular Value Decomposition) method can be used to obtain the deformation time series of the monitoring area, and the phase stacking (Stacking) technology can be used to obtain the deformation rate of the monitoring area.

[0145] The present invention also provides a specific embodiment, selecting the coastal area of ​​Xiamen as the experimental object, using 125 scenes of Sentinel 1A data, the experimental results are as follows: Figures 3 to 8 As shown, Figure 3 A time-space baseline diagram provided by an embodiment of the present invention, wherein the horizontal axis represents the date and the vertical axis represents the vertical baseline; Figure 4 A schematic diagram of the cumulative DEM error provided by an embodiment of the present invention, where different color depths represent different DEM errors; Figure 5 for Figure 4 Time series of DEM errors at points D1, D2, D3, and D4. The horizontal axis represents the date, the left vertical axis represents the accumulated vertical deformation, and the right vertical axis represents the DEM error. Figure 6 The DEM error deformation time series of the main regions provided by the embodiment of the present invention; Figure 7 for Figure 4 Comparison chart of DEM error correction and original time series at point D3; Figure 8 for Figure 4 Comparison chart of DEM error correction and original time series at point D4.

[0146] The InSAR time series DEM error estimation method provided by the present invention can not only obtain the deformation time series and deformation rate of the monitoring area, but also obtain the elevation time series of the time series ground landfill, and use the obtained elevation deformation time series to correct the DEM error in the deformation time series, thereby improving the accuracy of deformation monitoring.

[0147] The above are only a few embodiments of the present application and do not constitute any form of limitation to the present application. Although the present application is disclosed as above with preferred embodiments, it is not intended to limit the present application. Any technician familiar with the profession, without departing from the scope of the technical solution of the present application, using the technical content disclosed above to make slight changes or modifications are equivalent to equivalent implementation cases and fall within the scope of the technical solution.

Claims

1. An InSAR time series DEM error estimation method, characterized in that: The method comprises: S1, obtaining a plurality of single-view complex images of the area to be monitored, and performing differential interference processing and phase optimization processing on the plurality of single-view complex images to obtain a plurality of optimized interference patterns; S2, performing phase unwrapping on the multiple interference patterns to obtain multiple unwrapping patterns; S3, inverting the plurality of unwrapped images into a plurality of continuous phase images using an inversion function, and obtaining n grouped time points of the deformation phase time series of each pixel according to the plurality of continuous phase images; S4, using the n grouping time points to divide the deformation phase time series of the corresponding pixels into n+1 groups, and obtaining the DEM error corresponding to each group of pixels; The S2 specifically includes: S21, performing MCF phase unwrapping on the optimized multiple interference patterns to obtain multiple unwrapped patterns; S22, performing unwrapping error identification and phase correction on each pixel in the multiple unwrapping images to obtain multiple unwrapping images after phase correction; The S22 specifically includes: S221, normalizing corresponding pixels in a plurality of unwrapped images to obtain deformation phase rates of the pixels; S222, sorting the deformation phase rates of the pixels according to the middle time points of the corresponding interference patterns, and obtaining the median rate of the deformation phase rates of the sorted pixels in a time window of a preset time length before and after the current moment; S223, using the median rate to perform phase correction on the deformation phase rate of the pixel at the current moment, to obtain a plurality of phase-corrected unwrapping images.

2. The method according to claim 1, characterized in that: In S1, differential interference processing and phase optimization processing are performed on the plurality of single-view complex images to obtain a plurality of optimized interference patterns, which specifically includes: S11, preprocessing the plurality of single-view complex images, and performing differential interference processing on the preprocessed images in pairs to obtain a plurality of interference patterns; S12, using SRTM-DEM to remove the external terrain phase in the interference pattern to obtain a processed interference pattern; S13, extracting homogeneous points in the processed interference pattern, and using the homogeneous points to perform phase noise reduction on the interference pattern to obtain a plurality of optimized interference patterns.

3. The method according to claim 2, characterized in that The S3 specifically includes: S31, selecting an unwrapping graph that meets preset requirements from the multiple unwrapping graphs after phase correction; S32, using an inversion function to invert the unwrapped image that meets the preset requirements into a plurality of continuous phase images; S33, performing denoising on the multiple continuous phase images using a one-dimensional wavelet transform method to obtain a denoised deformation phase time series of each pixel; S34. Obtain n grouping time points according to the deformation phase time series of each pixel.

4. The method according to claim 3, characterized in that The S34 is specifically: The first-order derivative of the deformation phase time series of each pixel is taken to obtain n peak values ​​in the deformation phase time series of each pixel, and the time points corresponding to the peak values ​​are the grouping time points.

5. The method according to claim 3, characterized in that: The S34 is specifically: S341, determining a sliding window from the deformation phase time series of each pixel, and fitting the deformation phase within the sliding window to obtain a fitting function; S342, predicting the deformation phase of a future time period adjacent to the sliding window according to the fitting function to obtain a predicted phase time series; S343: Obtain the measured phase time series of the future time period, calculate the standard deviation or variance between the measured phase time series and the predicted phase time series, and use the peak value of the standard deviation or variance as the grouping time point.

6. The method according to claim 1, characterized in that The S4 specifically includes: S41, using the n grouping time points to divide the deformation phase time series of the corresponding pixels into n+1 groups; S42, selecting the unwrapped phases that meet the preset conditions in each group of pixels, constructing a DEM solution equation for the corresponding pixel group according to the unwrapped phases that meet the preset conditions, and obtaining the DEM error corresponding to each group of pixels according to the DEM solution equation.

7. The method according to claim 1 or 6, characterized in that: After S4, the method further includes: S5. De-noise the DEM errors corresponding to each group of pixels, and generate a ground elevation landfill time series based on the de-noised DEM errors.

8. The method according to claim 1 or 6, characterized in that: After S4, the method further includes: S6, denoising the DEM error corresponding to each group of pixels, and removing the corresponding denoised DEM error from the unwrapped phase of the unwrapped image to obtain an unwrapped image containing only the deformation phase; S7. Obtaining the deformation time series and deformation rate of the monitoring area according to the unwrapped image containing only the deformation phase.

Images