A PS-InSAR multi-baseline time-series phase unwrapping method based on a uniform-velocity deformation assumption

Abstract

This invention belongs to the field of synthetic aperture radar interferometry technology, specifically a PS-InSAR multi-baseline temporal phase unwrapping method based on the assumption of uniformly accelerated deformation. The method includes: preprocessing multi-temporal SAR data to obtain a temporal differential interferogram; selecting permanent scattering point (PS) on the temporal interferogram; constructing a second-order phase gradient constraint model in the time dimension under the assumption of uniformly accelerated deformation; using a two-way time-progression algorithm to solve the integer ambiguity gradient of the second-order phase gradient constraint model; and constructing and solving a globally weighted L1 norm optimization problem with respect to the absolute phase vector based on the integer ambiguity gradient. This invention, by establishing a second-order phase gradient constraint model, effectively describes the acceleration components in the surface deformation process, improving the adaptability to complex temporal deformation.

Description

Technical Field

[0001] This invention belongs to the field of synthetic aperture radar interferometry technology, and particularly relates to a PS-InSAR multi-baseline time-series phase unwrapping method based on the assumption of uniformly variable speed deformation. Background Technology

[0002] Traditional PS-InSAR time-series unwrapping methods typically construct single-baseline phase unwrapping (PU) models based on the Itoh condition assumption. These methods assume small phase changes between adjacent time steps, allowing for continuous reconstruction of the temporal phase through the accumulation of time-series phase differences. However, in regions with large-gradient surface deformation or nonlinear deformation processes, phase changes are drastic and often accompanied by frequent jumps, making it difficult for the continuity assumption to hold. This results in problems such as the Itoh condition not being met or error accumulation in the unwrapping results.

[0003] To improve unwrapping accuracy, previous studies have proposed temporal PU models based on multi-temporal baselines (MB). These models constrain the temporal phase evolution of the same pixel by combining multiple interferograms, thereby reducing the dependence of single-baseline methods on the Itoh condition. While the MB PU strategy improves PU accuracy to some extent, these methods typically employ the assumption of a constant deformation rate, meaning the deformation rate at the target point remains constant over time. This assumption is difficult to accurately characterize the actual deformation process in scenarios where the surface deformation rate varies significantly, exhibits nonlinear evolution, or undergoes abrupt changes, thus limiting the consistency between PU accuracy and temporal consistency. Summary of the Invention

[0004] The purpose of this invention is to address the problems existing in the prior art by proposing a PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly accelerated deformation assumption. This method introduces a second-order phase gradient constraint in the time dimension to characterize the acceleration term of the deformation phase over time, thereby modeling the nonlinear temporal deformation process. Based on this, this paper constructs a PS-InSAR multi-baseline temporal phase unwrapping model based on the uniformly accelerated deformation assumption and employs the minimum-cost network flow (MCF) method to globally optimize the temporal deformation phase. The proposed method can achieve robust unwrapping for large gradients and nonlinear deformations without relying on the Itoh condition.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A PS-InSAR multi-baseline temporal phase unwrapping method based on the assumption of uniformly accelerated deformation includes the following steps:

[0007] S1. Acquire multi-temporal SAR data of the same region and preprocess it to obtain a time-series differential interferogram;

[0008] S2. Select the permanent scatterer PS point on the time-difference interferogram.

[0009] S3. Construct a spatial triangulation based on the selected PS points to obtain the triangulation edge set; for each edge of the triangle, establish a deformation phase gradient model on each edge set based on the assumption of uniformly accelerated deformation; on this basis, introduce the second-order time phase gradient constraint to construct a second-order phase gradient constraint model.

[0010] S4. Using the main image time as a reference, a two-way time advancement algorithm is used to calculate the integer fuzziness gradient of the second-order phase gradient constraint model, and obtain the defuzzified temporal unwrapping phase gradient field.

[0011] S5: Using the defuzzified temporally unwrapped phase gradient field as the observation, construct and solve a global weighted L1 norm optimization problem with respect to the absolute phase vector, thereby obtaining the unwrapped phase of all PS points in the interferogram at all times.

[0012] By adopting the above technical solution, the present invention has the following beneficial effects:

[0013] 1. This invention introduces the assumption of uniformly accelerated deformation and establishes a second-order phase gradient constraint model in the time dimension, effectively describing the acceleration components in the surface deformation process, thereby accurately characterizing the nonlinear temporal deformation features. Compared with the traditional uniformly accelerated deformation assumption model, it significantly improves the adaptability to complex temporal deformation.

[0014] 2. By establishing a dual-objective minimization constraint, this invention simultaneously optimizes the integer fuzzy number gradient relationship of multiple time phases from a global perspective, significantly reducing the dependence on the Itoh condition and improving the temporal consistency of the unwrapping results.

[0015] 3. This invention effectively suppresses the uncertainty caused by phase abrupt changes through a bidirectional time-progression algorithm. Furthermore, by working in conjunction with the MCF method in the spatial dimension, it can achieve continuous and accurate unwrapping of large-gradient surface deformation regions, thus solving the problem of traditional temporal phase unwrapping failing in drastic deformation scenarios.

[0016] 4. This invention provides a PS-InSAR multi-baseline temporal phase unwrapping method based on the assumption of uniformly variable speed deformation. It can directly replace the temporal phase unwrapping step in the existing PS-InSAR processing flow, without the need for additional observation data or complex external constraints. It has good engineering feasibility and application value in nonlinear deformation monitoring scenarios. Attached Figure Description

[0017] Figure 1 This is the overall flowchart of the present invention;

[0018] Figure 2 This is a schematic diagram of the data preparation stage of the present invention; where (a) is the reference absolute phase of the simulation, (b) is the winding phase of the simulation, (c) is the scatter plot of the expansion center of the simulation as a function of volume, and (d) is the line graph of the triangular mesh edge replacement effect as a function of K value.

[0019] Figure 3 This is a schematic diagram of the triangular mesh constructed by the selected PS in this invention; where (a) is the Delaunay triangular mesh, and (b) is the triangular mesh after filtering and retaining the highly coherent edges of the Delaunay triangular mesh and using KNN to fill the edges. The K value is set to 100 in the figure.

[0020] Figure 4 The diagrams show the time-series unwrapping errors using different methods; where (a) represents the time-series unwrapping error of the Itoh+MCF+Delaunay method, (b) represents the time-series unwrapping error of the Itoh+MCF+Delaunay, KNN method, (c) represents the time-series unwrapping error of the MB PU (uniform velocity assumption)+MCF+Delaunay method, (d) represents the time-series unwrapping error of the MB PU (uniform velocity assumption)+MCF+Delaunay, KNN method, (e) represents the time-series unwrapping error of the MB PU (uniform velocity assumption)+MCF+Delaunay method, and (f) represents the time-series unwrapping error of the MB PU (uniform velocity assumption)+MCF+Delaunay, KNN method.

[0021] Figure 5Line plots showing the temporal phase fitting of the PS at the center of the scene for different methods are provided. Among them, (a) is the temporal phase fitting line plot of the Itoh+MCF+Delaunay method, (b) is the temporal phase fitting line plot of the Itoh+MCF+Delaunay, KNN method, (c) is the temporal phase fitting line plot of the MB PU (uniform velocity assumption)+MCF+Delaunay method, (d) is the temporal phase fitting line plot of the MB PU (uniform velocity assumption)+MCF+Delaunay, KNN method, (e) is the temporal phase fitting line plot of the MB PU (uniform velocity assumption)+MCF+Delaunay method, and (f) is the temporal phase fitting line plot of the MB PU (uniform velocity assumption)+MCF+Delaunay, KNN method. Detailed Implementation

[0022] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings.

[0023] like Figure 1 As shown in the figure, this embodiment provides a PS-InSAR multi-baseline temporal phase unwrapping method based on the assumption of uniformly variable speed deformation, which includes the following steps:

[0024] S1. Acquisition of multi-temporal SAR data: Acquire multi-temporal SAR data from the same region for temporal phase unwrapping.

[0025] S2. Preprocessing of multi-temporal SAR data. This step specifically includes:

[0026] S21. From the multi-temporal SAR data, select one image as the common master image, and define the other images in the multi-temporal SAR data other than the common master image as slave images. In this embodiment, SNAP software is used to select the common master image, and the image located at the middle date of the time series is used as the master image.

[0027] S22. Register each secondary image with the common primary image to generate a series of registered primary-secondary image pairs.

[0028] S23. Perform differential interferometry on each master-slave image pair to generate an initial differential interferogram; and sequentially remove the flat phase, orbital error phase, and terrain phase from the initial differential interferogram to obtain a temporal differential interferogram.

[0029] S3. Construct a second-order phase gradient constraint model by introducing a second-order phase gradient in the time dimension for constraint; this step specifically includes:

[0030] S31. Construct a spatial triangulation based on the selected PS points and remove duplicate edges to obtain the initial triangulation edge set. To ensure the connectivity of the triangulation, this embodiment further optimizes the edge set after obtaining it. The optimization process is as follows:

[0031] S311. Based on the time series phase data of the two PS points corresponding to each edge, calculate the temporal coherence between the PS point pair along all interferograms;

[0032] S312. Remove edges with temporal coherence below a preset threshold from the edge set;

[0033] S313. For isolated PS points or subnetworks caused by edge removal, the K-Nearest Neighbor (KNN) algorithm is used to find new connecting edges, and the new edges that meet the time coherence threshold requirement are added to the edge set to form the final triangular network edge set.

[0034] S32. Based on the optimized triangular mesh edge set from step S31, for each edge of the triangle, based on the assumption of uniformly accelerated deformation, establish a deformation phase gradient model on each edge set; on this basis, introduce a second-order time phase gradient constraint to construct a second-order phase gradient constraint model. Specifically:

[0035] S321. Assuming the phase corresponding to the deformation is generated by a uniformly accelerated process, in the... The absolute phase (deformation phase) at each PS point satisfies the following formula

[0036] (1);

[0037] in This represents the difference between the starting image time and the main image time, which will be expressed as follows: Since all images are registered with the common master image, the initial time will be used here. Set to the time corresponding to the main image. for Time of the first Absolute phase at each PS point, Main image moment Absolute phase at each PS point, Indicates the radar wavelength. Indicates the main image time. The initial deformation velocity at point PS express Time of the first Deformation acceleration at point PS.

[0038] S322. For each edge in the triangular mesh edge set, based on the uniformly accelerated deformation assumption of formula (1), the deformation phase model of point PS is converted into a deformation phase gradient model on the edge. The expression of this deformation phase gradient model is:

[0039] (2);

[0040] in, Given an edge on the set of edges of a triangulation network, express The phase gradient of deformation at any moment, This indicates the difference in deformation rate at different times in the main image. express Difference in deformation acceleration at any given moment.

[0041] S323. Introduce a second-order phase gradient constraint in the time dimension and construct a second-order phase gradient constraint model; specifically:

[0042] 1) Define the mathematical relationship between the wrapped phase and the absolute phase to be solved as follows:

[0043] (3);

[0044] in, for Time of the first Absolute phase at each PS point, for Time of the first The entanglement phase at each PS point for Time of the first The number of integer blurs at each PS point;

[0045] 2) Based on the mathematical relationship between the wound phase and the absolute phase to be solved, the deformation phase gradient model is rewritten as the wound phase gradient plus... Product with the gradient of the integer fuzzy number:

[0046] (4);

[0047] in, express The phase gradient is constantly wrapped around. express Gradient of integer fuzzy number at time intervals;

[0048] 3) By rearranging the above equation and taking the deformation acceleration difference and integer ambiguity gradient corresponding to the triangular mesh edges at a given time as the unknowns in the equation, a second-order phase gradient constraint model is obtained:

[0049] (5).

[0050] S4. Using the main image time as a reference, a two-way time-advancement algorithm is employed to calculate the integer fuzziness gradient of the second-order phase gradient constraint model, thereby obtaining the defuzzified temporal unwrapping phase gradient field. This step specifically includes:

[0051] S41. Using the main image time as the zero constraint center, establish a bi-objective minimization constraint to constrain the continuity of the deformation acceleration difference corresponding to the triangular mesh edge set:

[0052] (6)

[0053] ;

[0054] in, Indicates the moment of the main image. For the set of edges of the triangular mesh, The search range for integers is typically set to [-1, 1];

[0055] S42, with the main image time As a zero-phase-difference reference moment, the initial velocity difference at this moment is set. The value is 0; for the deformation acceleration difference corresponding to the triangular mesh edge set, a bidirectional time-advancement algorithm is used to advance from the main image time to both the preceding and following time steps of the time series. Combined with the dual-objective minimization constraint in step S41, the second-order phase gradient constraint model in step S3 is solved to obtain the integer ambiguity gradient candidate values ​​of the interferogram at each time step; in each advancement direction, for the current time step, the candidate values ​​of the integer ambiguity gradient are searched within the preset integer candidate set. Compared with the unidirectional advancement method starting from the first time step, the bidirectional time-advancement algorithm used in this embodiment divides the original full-length T cumulative error path into two segments, each with a length of approximately T / 2. Theoretically, this can reduce the variance of the time cumulative error by about 1 / 4, ensuring the continuity of the constraint at the main image, eliminating boundary effects, and avoiding the "blind initialization" problem of the unidirectional algorithm at the first time step, i.e., the first time step t1 lacks the previous time step t0 used for difference, resulting in a lack of time constraints.

[0056] S43. Based on the dual-objective minimization constraint obtained in step S41, select the candidate value that satisfies the dual-objective minimization constraint from the candidate values ​​of the integer fuzzy number gradient obtained in step S42 as the optimal integer fuzzy number gradient; substitute this optimal integer fuzzy number gradient into the relationship between the absolute phase gradient and the entanglement phase gradient to solve for the estimated unentanglement phase gradient. .

[0057] S5. Using the deblurred temporally unwrapped phase gradient field as the observation, construct and solve a globally weighted L1 norm optimization problem with respect to the absolute phase vector, thereby obtaining the unwrapped phase of all PS points in the interferogram at all times. This step specifically includes:

[0058] S51, Select A permanent scatterer, based on A triangulation is constructed using a permanent scatterer, and the triangulation is the first... An edge is represented by its start point and end point as follows: S52. The global weighted L1 norm optimization problem is expressed as:

[0059] ;

[0060] ;

[0061] ;

[0062] ;

[0063] in, It is the temporal coherence weighting coefficient. It is the total number of edges of the triangular network. X is the gradient operation matrix; X is the absolute phase vector to be solved, where each element represents... b is the absolute phase to be determined for a single PS point at time t; b is the vector corresponding to the deblurred temporal unwrapped phase gradient field, where each element represents The unwrapped phase gradient on the edge of a single triangular mesh at any given time; Indicates the first The triangulation in the interferogram The winding phase gradient on the edge, This indicates the total number of interferograms.

[0064] S53. Construct the gradient operation matrix A using PS point pairs on the edges of the triangular mesh, and solve the global weighted L1 norm optimization problem with respect to the absolute phase vector to obtain the final temporal phase unwrapping result. Wherein:

[0065] The matrix A is:

[0066] ;

[0067] Among them, matrix Each row corresponds to an edge on the triangular network, and the matrix... Each column corresponds to a PS point, for the th in the triangulation For each row corresponding to an edge: the PS position at the edge's starting point has a value of 1, the PS position at the edge's ending point has a value of -1, and all other positions have a value of 0; the algorithm complexity of this solution process is O(n log n). The overall algorithm complexity is .

[0068] Experimental verification:

[0069] like Figure 2 As shown, the volume change at the expansion center in a uniform elastic half-space at a depth of 2 km was simulated, and 4000 sampling points were randomly selected as PS points within a 20 km × 20 km area. Assuming a radar wavelength of 0.056 m and a satellite repetition period of 35 days, the relative deformation occurring at these PS points was calculated using 20 randomly spaced satellite passes, and converted into the line-of-sight phase difference relative to the 10th time step, thus obtaining the reference absolute phase of the time series. Finally, using the 10th time step as the master image time, phase difference calculations were performed on the remaining slave images and then wrapped with the master image to obtain the wrapped phase of the time series. The simulated absolute phase is shown below. Figure 2 As shown in (a), the winding phase is as follows Figure 2 As shown in (b), the volume change at the expansion center is as follows: Figure 2 As shown in (c), the deformation mainly occurs at the center of the scene. At time 10, the difference between the main image and itself is calculated, resulting in an absolute phase of zero. Since we are concerned with the relative deformation between the image time and the main image time, the absolute phase at the main image time is not considered in the subsequent recording of the unwrapping results. In the experiment, to facilitate the evaluation of the error between the unwrapping result and the reference absolute phase, the phase value corresponding to point PS at the upper right corner of the scene is used as a reference. The difference between the unwrapped phase and the reference absolute phase at point PS is calculated, and this residual value is subtracted from the unwrapped phase to obtain the final unwrapped phase for result comparison.

[0070] The unwrapping process of temporal differential interferograms (TIAs) consists of two steps: one-dimensional (1D) temporal unwrapping and two-dimensional (2D) spatial sparse grid unwrapping. In the 1D unwrapping step, the methods compared are: constructing a phase gradient field based on the Itoh condition; using a multi-baseline unwrapping method based on the uniform deformation assumption (MBPU, uniform deformation assumption) to estimate the unwrapped phase gradient; and using the multi-baseline unwrapping method based on the uniformly variable deformation assumption (MBPU, uniformly variable deformation assumption) employed in this invention to estimate the unwrapped phase gradient. In the 2D unwrapping step, the MCF method based on the L1 norm is used to unwrap the differential interferograms at each time step. During the triangulation construction process, depending on whether the Delaunay triangulation edge set is optimized, the methods compared employ the traditional Delaunay triangulation and a method using a KNN network to optimize the triangulation construction, abbreviated as Delaunay, KNN. The effect of triangulation edge replacement varies with the value of K in the KNN as shown below. Figure 2 As shown in Figure (d), the change in average coherence after edge replacement is statistically analyzed with a step size of 50. When K=100, the coherence increases the fastest. After K>100, the increase in coherence slows down. The larger the value of K, the more edges are retained and the longer the computation time. In order to balance computational efficiency and performance gain, the value of K is set to 100 in the experiment.

[0071] like Figure 3 As shown, a triangular mesh is constructed based on the selected PS point, where the solid black circles represent the selected PS points and the short black lines represent the edges of the triangular mesh. Figure 3 In (a), the number of edges in the Delaunay triangulation constructed based on the selected PS points after removing duplicate edges is 11974. Figure 3 The number of sides in the triangular mesh after edge screening and patching in (b) is 11568.

[0072] To measure the degree of approximation between the unwrapping result and the reference absolute phase, the root mean square error (RMSE) was chosen as the quantification metric. The lower the RMSE, the higher the unwrapping accuracy. Figure 4 In (a), the method constructs a phase gradient field using the Itoh condition and performs MCF unwrapping at each time step. However, significant unwrapping errors occur at times with large deformation gradients, i.e., when significant deformation occurs. This is because the method is limited by the Itoh condition, leading to decreased unwrapping accuracy in scenarios with phase discontinuities. Furthermore, this method does not consider the dependency between images from adjacent time steps, resulting in larger unwrapping errors at these times. Figure 5 The trend of the line chart of the untangling phase in (a) is inconsistent with that of the line chart of the reference absolute phase. Figure 4 (b) in Figure 4 Based on method (a), this method filters and retains edges with high coherence in the Delaunay triangulation and uses KNN to complete the unconnected PS, significantly improving the tolerance of this single-baseline method in phase discontinuity regions and achieving accurate unwrapping in most cases. However, when the deformation gradient is large and the Itoh condition is not satisfied, this method still cannot accurately unwrap the network or capture the data. Figure 5 The trend of change in the reference absolute phase in the piecewise linear plot at the moment of large gradient deformation (b). Since the MB PU method does not depend on the Itoh condition, therefore... Figure 5 The method based on the uniform velocity assumption in (c) can generally capture the temporal variation trend of the reference absolute phase, but due to the uniform deformation assumption, there are still significant errors at large gradient deformation moments. After coherence screening of the triangular mesh edge set, the temporal piecewise linear fitting of this method is further improved, but there are still some differences from the temporal variation of the reference absolute phase. Figure 4 In (e), after using the MB PU method based on the uniformly variable speed assumption proposed in this invention, the unwrapping error is significantly lower than that of other MB PU unwrapping strategies. Figure 4 (c) and Figure 4 A significant decrease was achieved in (d), in Figure 5 Visually, except for times 3, 4, and 5, the broken lines at other times largely coincide with the reference absolute phase broken line, and the error is relatively small at times 3, 4, and 5. After coherence screening of the triangular mesh edge set, the RMSE value of the unwrapping result of the method proposed in this invention is near 0 at all times, and the trend of the broken line graph of the unwrapped phase highly coincides with the reference absolute phase, proving that the unwrapping performance of the proposed method in the large gradient nonlinear deformation region is better than the comparative method.

[0073] Table 1: Goodness of fit (R²) of temporal phase-reference absolute phase at scene center PS for different unwrapping methods

[0074]

[0075] This embodiment quantifies the effectiveness of different unwrapping methods, using the coefficient of determination (R²) as the evaluation index to analyze the degree of fit between the temporal phase at the PS position of the scene center and the reference absolute phase. The R² results for each method are listed in Table 1. R² values ​​are less than or equal to 1; the closer R² is to 1, the higher the degree of fit. It should be noted that since the reference absolute phase in the simulation experiment does not include noise, terrain phase, or atmospheric phase, theoretically, the degree of fit between the unwrapping result and the reference absolute phase can reach 1. As can be seen from Table 1, without performing a screening operation on the triangular mesh edge set, the R² of the single-baseline unwrapping method is much lower than that of the multi-baseline unwrapping method based on the uniform velocity assumption and the uniformly variable velocity assumption. This is due to the limitation of the Itoh condition, which causes a decrease in unwrapping accuracy in the large gradient deformation region, resulting in a lower degree of fit of the temporal piecewise linear plot. Furthermore, the MB PU (uniformly variable velocity assumption) + MCF method proposed in this invention can achieve an R² higher than 0.99 even without performing a coherence screening operation on the triangular mesh edge set, significantly outperforming the comparative methods. After filtering and retaining highly coherent edges and performing edge-filling operations using KNN, only the method proposed in this invention achieved R² of 1. This indicates that the MB PU method based on the uniformly variable speed assumption estimates the unwrapped phase gradient with higher accuracy, and the unwrapped phase estimated based on this gradient field has the best fit with the reference absolute phase. It also has high unwrapping accuracy in the large gradient nonlinear deformation region, and its unwrapping performance is significantly better than the comparative methods.

[0076] Of course, the present invention may have other various embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and modifications according to the present invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A PS-InSAR multi-baseline temporal phase unwrapping method based on the assumption of uniformly variable speed deformation, characterized in that, Includes the following steps: S1. Acquire multi-temporal SAR data of the same region and preprocess it to obtain a time-series differential interferogram; S2. Select the permanent scatterer PS point on the time-difference interferogram. S3. Construct a spatial triangulation network based on the selected PS points to obtain the set of edges of the triangulation network; for each edge of the triangle, based on the assumption of uniformly accelerated deformation, establish a deformation phase gradient model on each edge set; on this basis, introduce a second-order time phase gradient constraint to construct a second-order phase gradient constraint model; the specific steps for constructing the second-order phase gradient constraint model include: S31. Convert the deformation phase model at point PS into a deformation phase gradient model on the edge, the expression of which is: ； in, Given an edge on the set of edges of a triangulation network, express The phase gradient of deformation at any moment, Indicates the radar wavelength. This represents the difference in deformation rate at the initial moment. Indicates the current moment. express Difference in deformation acceleration at any moment; S32. Introduce a second-order phase gradient constraint in the time dimension to construct a second-order phase gradient constraint model, the expression of which is: ； in, express The phase gradient is constantly wrapped around. express Gradient of fuzzy number at time intervals; S4. Using the main image time as a reference, a two-way time advancement algorithm is used to calculate the integer fuzziness gradient of the second-order phase gradient constraint model, and obtain the defuzzified temporal unwrapping phase gradient field. S5. Using the defuzzified temporally unwrapped phase gradient field as the observation, construct and solve a global weighted L1 norm optimization problem about the absolute phase vector, thereby obtaining the unwrapped phase of all PS points in the interferogram at all times.

2. The PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly variable deformation assumption as described in claim 1, characterized in that, The preprocessing process in step S1 includes: S11. Select one image from the multi-temporal SAR data as the common master image, and define the other images in the multi-temporal SAR data other than the common master image as slave images. S12. Register each scene image with the common master image to generate a series of registered master-slave image pairs. S13. Perform differential interferometry on each master-slave image pair to generate an initial differential interferogram; and sequentially remove the flat phase, orbital error phase, and terrain phase from the initial differential interferogram to obtain a temporal differential interferogram.

3. The PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly variable deformation assumption as described in claim 1, characterized in that... The process of constructing a spatial triangulation based on the PS point and obtaining the triangulation edge set in step S3 includes the step of removing duplicate connected edges after constructing the spatial triangulation.

4. The PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly variable speed deformation assumption as described in claim 1, characterized in that... Step S3 further includes optimizing the triangulation edge set after obtaining it: Based on the time series phase data of the two PS points corresponding to each edge, calculate the temporal coherence between the PS point pair along all interferograms; Edges with temporal coherence below a preset threshold are removed from the edge set; For isolated PS points or subnetworks resulting from edge removal, the K-Nearest Neighbor (KNN) algorithm is used to find new connecting edges, and new edges that meet the time coherence threshold requirement are added to the edge set to form the final triangulation edge set.

5. The PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly variable deformation assumption as described in claim 1, characterized in that... The process of introducing a second-order phase gradient constraint in the time dimension to construct a second-order phase gradient constraint model in step S32 also includes a model arrangement step that treats the deformation acceleration difference as an unknown quantity, specifically including the following sub-steps: S321. The mathematical relationship between the entangled phase and the absolute phase to be solved is defined as follows: ； in, for Time of the first Absolute phase at each PS point, for Time of the first The entanglement phase at each PS point for Time of the first The number of integer blurs at each PS point; S322. Based on the mathematical relationship between the wrapped phase and the absolute phase to be solved, the deformation phase gradient model is rewritten as the wrapped phase gradient plus... Product with the gradient of the integer fuzzy number: 。 6. The PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly variable deformation assumption as described in claim 5, characterized in that... Step S4, based on the main image time, employs a bidirectional time-progression algorithm to calculate the integer fuzziness gradient of the second-order phase gradient constraint model, thereby obtaining the defuzzified temporal unwrapping phase gradient field; specifically: S41. Using the main image time as the zero constraint center, establish a bi-objective minimization constraint to constrain the continuity of the deformation acceleration difference corresponding to the triangular mesh edge set: in, Indicates the moment of the main image; The set of edges of a triangular mesh; S42, with the main image time As a zero-phase-difference reference moment, the initial deformation rate difference at this moment is set. =0; For the deformation acceleration difference corresponding to the edge set of the triangular network, starting from the main image time, it is simultaneously advanced to the preceding and following time of the time series. Combined with the dual objective minimization constraint of step S41, the second-order phase gradient constraint model of step S3 is solved to obtain the integer ambiguity gradient candidate values ​​of the interferogram at each time. In each advancing direction, for the current time, the candidate values ​​of the integer ambiguity gradient are searched in the preset integer candidate set [-1,1]. S43. Based on the dual-objective minimization constraint obtained in step S41, select the candidate value that satisfies the dual-objective minimization constraint from the candidate values ​​of the integer fuzzy number gradient obtained in step S42 as the optimal integer fuzzy number gradient; substitute this optimal integer fuzzy number gradient into the relationship between the absolute phase gradient and the entanglement phase gradient to solve for the estimated unentanglement phase gradient. .

7. A PS-InSAR multi-baseline temporal phase unwrapping method based on the uniformly variable deformation assumption as described in claim 6, characterized in that... The specific steps of step S5 include: S51, Select A permanent scatterer, based on A triangulation is constructed using a permanent scatterer, and the triangulation is the first... An edge is represented by its start point and end point as follows: ; S52. The global weighted L1 norm optimization problem can be expressed as: ； ； ； ； in, It represents the total number of edges in the triangular mesh; X is the gradient operation matrix; X is the absolute phase vector to be solved, where each element represents... b is the absolute phase to be determined for a single PS point at time t; b is the vector corresponding to the deblurred temporal unwrapped phase gradient field, where each element represents The unwrapped phase gradient on the edge of a single triangular mesh at any given time; It is the temporal coherence weighting coefficient; Indicates the first The triangulation in the interferogram The winding phase gradient on the edge, Indicates the total number of interferograms; S52. Construct the gradient operation matrix using PS point pairs on the edges of the triangular mesh. Solve the globally weighted L1 norm optimization problem with respect to the absolute phase vector; where: The matrix for: ； Among them, matrix Each row corresponds to an edge on the triangular network, and the matrix... Each column corresponds to a PS point, for the th in the triangulation The row corresponding to each edge: the PS position at the start of the edge has a value of 1, the PS position at the end of the edge has a value of -1, and the PS position at other positions has a value of 0.

Images