Ground surface deformation monitoring method based on intermittent stacking time sequence InSAR

Abstract

The invention belongs to the technical field of terrain remote sensing monitoring, and particularly relates to a surface deformation monitoring method based on an intermittent stacking time sequence InSAR (Intermittent Synthetic Aperture Radar), which comprises the following steps of: firstly, acquiring a plurality of SAR images, high-precision terrain data and orbit data of a target area in a lifting orbit within a period of time; performing geographical registration, interference combination selection, space-time filtering and phase unwrapping on the SAR image to generate a time sequence interferogram; calculating an annual average deformation rate of each pixel point in the interferogram, and generating an annual average deformation rate result of the target area; and finally, generating a suspicious landslide polygon according to an annual average deformation rate result, and performing visual interpretation and field investigation to obtain a real area of the landslide polygon. According to the method, result missing caused by the intermittent coherence characteristic of the pixel points can be avoided, the pixel selection has higher latitude, higher coverage can be achieved, the requirement for computing power is remarkably reduced, and a deformation rate result in a low-coherence environment can be rapidly obtained.

Description

Technical Field

[0001] The invention belongs to the technical field of terrain remote sensing monitoring, and in particular relates to a surface deformation monitoring method based on intermittent stacking time series InSAR. Background Art

[0002] InSAR (Interferometric Synthetic Aperture Radar) technology is an active microwave remote sensing technology with outstanding advantages such as non-contact, all-day, all-weather, high precision and high resolution. It can obtain deformation data on a large range of the surface and is used in active landslide detection, pre-landslide prediction, landslide range and volume estimation, landslide risk assessment and other fields. It has become the most widely used technology in the field of landslide monitoring today.

[0003] Time-series InSAR technology can effectively monitor surface deformation. However, due to the steep terrain and dense vegetation coverage in river canyon areas, the coherence of short-wave SAR satellites such as Sentinel-1 is very low in these areas. Especially in the long-time series landslide monitoring scenario, the acquired images have intermittent coherence (intermittent coherence refers to the property that ground scatterers are difficult to maintain high coherence in a long time series). As a result, the traditional time-series InSAR method can obtain very few deformation results in low-coherence areas, making it impossible to conduct effective landslide monitoring and identification.

[0004] Stacking-InSAR (Stacked Synthetic Aperture Radar Interferometry) and SBAS-InSAR (Small Baseline Subset Synthetic Aperture Radar Interferometry) technologies are the most commonly used time-series InSAR technologies in the field of geological disaster deformation monitoring. In 1998, Sandwell et al. proposed the Stacking-InSAR algorithm, which combines and weights M unwrapped interferometric phase images to weaken the influence of orbit, atmosphere and terrain errors in the interferometric phase. The specific formula for calculating the annual average deformation rate is as follows:

[0005]

[0006] Among them, V LOS is the annual average deformation rate calculated by Stacking, λ is the radar wavelength, ΔΦ is the phase difference vector, ΔT is the time span vector, <ΔΦ, ΔT> is the vector dot product, M is the total number of interference pairs, is the time series interference diagram after unwrapping; ΔT M is the time span of the interference pattern.

[0007] The basic principle of SBAS-InSAR technology is to minimize the errors and spatiotemporal incoherence caused by the difference in viewing angles between interference pairs by setting the spatiotemporal baseline threshold, while at the same time using high-pass filtering in the time domain and low-pass filtering in the space domain to weaken the influence of the atmospheric phase.

[0008] However, in the application scenarios of landslide monitoring and identification in long time series and low coherence areas, few pixels can maintain high coherence in long time series. As long as there is a missing pixel in any interference pattern, the pixel will be missing in the result of the entire time series. This is a defect of traditional SBAS-InSAR and Stacking-InSAR time series algorithms.

[0009] In 2007, Biggs et al. first introduced the intermittent coherence method to improve the estimation of fault deformation linear velocity. In 2013, Sowter et al. introduced the intermittent coherence method based on the SBAS algorithm, increasing the solution scale of the SBAS algorithm from the entire image to the pixel level. The proposed ISBAS algorithm can obtain a result coverage range 4-40 times that of the SBAS algorithm. However, in the actual deformation monitoring and geological disaster monitoring needs, the SBAS algorithm has a long calculation time and high requirements for computer running memory, and the ISBAS algorithm consumes a lot of time and computing resources.

[0010] An existing patent with patent application number CN202310054817.6 discloses an InSAR time series deformation monitoring method with automatic error correction. The scheme discloses stacking time series deformation solution for a set of disentangled interference patterns after quality screening in an intermittent manner. However, the scheme screens the entire interference image set according to certain standards, and then artificially eliminates some interference pair images, which will waste potentially useful interference information.

[0011] An existing patent application with patent application number CN202310605697.4 discloses a method for correcting and screening time-series InSAR interseismic deformation interference patterns. The scheme uses an intermittent coherence method to perform phase estimation on interference pixels based on common superposition to obtain the disturbance error of each time phase. However, the scheme is mainly for a certain moment in the earthquake time, and selects SAR images for a period of time before and after the earthquake for interference, and uses the intermittent Stacking method to optimize the density of coherent points. The scheme mainly processes several interference pairs at a certain moment, not interference pairs spanning a long time series.

[0012] In view of the above defects of the existing technology, it is necessary to develop a terrain deformation monitoring method suitable for time-series InSAR in a low-coherence environment, which can meet the needs of landslide deformation monitoring in a long-term series and low-coherence environment. Summary of the invention

[0013] The present invention provides a surface deformation monitoring method based on intermittent stacking time series InSAR. The method selects a corresponding interference pattern for each pixel point, which can avoid the missing results caused by the intermittent coherence characteristics of the pixel points. In terms of pixel selection, missing values ​​are allowed in the interference patterns at some time points, and only those interference patterns with phase values ​​are considered. Compared with the existing ISBAS technology, the method has a higher tolerance and can achieve a higher coverage. The method can flexibly handle missing values ​​and low-coherence areas, significantly reduce the demand for computing power, and can achieve rapid acquisition of InSAR deformation rate results in a low-coherence environment, and achieve a rapid survey of landslide disasters in a wide-area low-coherence environment. It has great application value in the field of landslide and other geological disaster monitoring.

[0014] The method includes the following steps: data acquisition: obtaining multiple SAR images of the target area in ascending and descending orbits within a period of time, obtaining high-precision terrain data of the target area, and obtaining SAR orbit data. Standard differential interferometry processing: geo-registering the acquired SAR images, calculating the spatiotemporal baseline between each SAR image according to the SAR orbit data, setting the spatiotemporal baseline threshold according to the spatiotemporal baseline, selecting the interference combination based on the spatiotemporal baseline threshold, and generating a time series interference pair; performing spatiotemporal filtering and phase unwrapping on the generated time series interference pair to generate an unwrapped time series interference map. Intermittent stacking time series InSAR processing: performing unwrapping error and atmospheric delay checks on a certain pixel point in the unwrapped time series interference map, eliminating the interference map with the missing pixel point, and then performing intermittent stacking time series InSAR processing on the pixel point in the remaining interference map, and calculating the annual average deformation rate of the pixel point; repeating this step, calculating the annual average deformation rate of each pixel point in the time series interference map, and finally generating the annual average deformation rate result of the target area. Landslide hazard generation and cataloging: Based on the annual average deformation rate results of the target area, the deformation area is extracted and the suspected landslide polygons are generated; the generated suspected landslide polygons are visually interpreted and field investigated to obtain the true area of ​​the landslide polygon and catalog the landslide.

[0015] In a specific implementation, the GMTSAR open source software installed on the Ubuntu platform is used to perform the standard differential interferometry processing step, and the SBAS technology is used to select the interference combination in the standard differential interferometry processing step.

[0016] In a specific implementation, in the standard differential interferometry processing step, a Goldstein filter is used to perform spatiotemporal filtering on the generated time series interferometry pairs.

[0017] In a specific implementation, in the standard differential interference processing step, the average coherence of the water area is collected, and the average coherence of the water area is set as the threshold of interference phase unwrapping, and the interference phase is phase unwrapped to generate an unwrapped time series interference pair.

[0018] In a specific embodiment, in a standard differential interferometry processing step, the interferometric phase is phase unwrapped using Snaphu software.

[0019] In a specific embodiment, in the intermittent stacking time series InSAR processing step, the following calculation formula is used to calculate the annual average deformation rate V of each pixel point in the interference pattern: LOS (x, y):

[0020]

[0021] Wherein, (x, y) represents the pixel point with row and column coordinates x and y in the geographic coordinate system, ΔΦ(x, y) represents the phase difference vector; m(x, y) is the pixel point (x, y) in the mth interference pattern in the interference pair combination; is the phase value of the pixel point (x, y) in the mth interference pattern in the interference pair combination; m(x ref ,y ref ) is the reference point (x) in the mth interference pattern of the interference pair combination. ref ,y ref ); is the reference point (x ref ,y ref ) is the phase value in the mth interference pattern in the interference pair combination; ΔT(x, y) represents the time span vector of the pixel point (x, y); ΔT m(x，y) is the time span of the pixel point (x, y) in the mth interference pattern in the interference pair combination; ∑ represents a diagonal matrix; V LOS (x, y) is the annual average deformation rate corresponding to the pixel point (x, y); <ΔΦ(x, y), ΔT(x, y)> is the vector dot product.

[0022] In a specific implementation, in the intermittent stacking time series InSAR processing step, the number of interference pair combinations and the coefficient matrix rank of each pixel point are also calculated:

[0023] Noi(x,y)=m(x,y) (3)

[0024] Roi(x,y)=rank(B(x,y)) (4)

[0025] Wherein, Noi(x, y) is the number of interference pair combinations of pixel point (x, y), and Roi(x, y) is the coefficient matrix rank of pixel point (x, y).

[0026] In a specific implementation, in the step of generating and cataloguing landslide hazards, the deformation area is extracted by using mathematical statistics on the annual average deformation rate results, and the suspicious landslide polygons are identified and generated by using a spatial clustering method.

[0027] In a specific implementation, in the step of generating and cataloguing landslide hazards, the step of using mathematical statistics to extract deformation areas from the annual average deformation rate results is as follows: calculating the mean μ and standard deviation σ of the annual average deformation rate results, and based on the principle that the error part of the unwrapped phase diagram obeys the normal distribution, dividing the results outside the range of μ±3σ into deformation areas, and dividing the results within the range of μ±3σ into non-deformation areas.

[0028] In a specific implementation, in the step of generating and cataloguing landslide hazards, the step of using the spatial clustering method to identify and generate suspicious landslide polygons for the annual average deformation rate results is as follows: the first step of spatial clustering analysis is performed based on the hot spot analysis algorithm of the Getis-Ord Gi* statistic to obtain highly clustered feature points, and the calculation formula is as follows:

[0029]

[0030] Among them, n represents the number of all feature points, n ij represents the number of all feature points within the search distance d of the i-th feature point (i = 1, 2, ..., n; j = 1, 2, ..., n ij ; j≠i); x represents the value of the feature point; and S* represent the mean and standard deviation of all feature point values; the aggregation point calculation method is used to eliminate discrete points and extract aggregation surfaces. The boundaries of the aggregation surfaces are connected to obtain suspicious landslide polygons.

[0031] The present invention has at least the following beneficial effects:

[0032] 1. The method of the present invention can avoid the missing results caused by the intermittent coherence characteristics of the pixel points by selecting the corresponding interference pattern for each pixel point. In terms of pixel selection, missing values ​​are allowed in the interference patterns at some time points, and only those interference patterns with phase values ​​are considered. Compared with the existing ISBAS technology, it has higher tolerance and can achieve higher coverage.

[0033] 2. The method of the present invention significantly reduces the demand for computing power by flexibly processing missing values ​​and low-coherence areas, and can achieve rapid acquisition of InSAR deformation rate results in low-coherence environments.

[0034] 3. The method of the present invention can be used as a complete technical process method suitable for engineering practice, which can realize the rapid survey of landslide disasters in a wide-area low-coherence environment and has great application value in the field of geological disaster monitoring such as landslides. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is a flow chart of an embodiment of the present invention.

[0036] Figure 2 Schematic diagram of the principle of intermittent stacking time series InSAR processing.

[0037] Figure 3 Schematic diagram of highly clustered feature points obtained by hotspot analysis in an embodiment of the present invention.

[0038] Figure 4 It is a schematic diagram of the process from extracting to aggregating surfaces to generating suspected landslide polygons in an embodiment of the present invention. DETAILED DESCRIPTION

[0039] See also Figure 1 The present invention provides a surface deformation monitoring method based on intermittent stacking time series InSAR, which specifically includes the following steps:

[0040] Data acquisition: Sentinel SAR (synthetic aperture radar of the European Space Agency's Sentinel-1 satellite) is used to obtain multiple SAR images of the target area in ascending and descending orbits over a period of time, obtain high-precision terrain data of the target area, and obtain the orbital data of the Sentinel SAR, which includes information such as the spatial position and operating speed of the Sentinel SAR.

[0041] By collecting multiple SAR images in both ascending and descending directions, more comprehensive surface deformation information can be obtained. Images in different orbital directions can provide complementary information, which helps to improve the quality of the interferogram and the accuracy of the interpretation. High-precision terrain data refers to the digital elevation model (DEM), which provides detailed terrain information in the study area, including altitude, slope, etc. Since the propagation path of SAR signals will be affected by the undulating terrain, an accurate DEM is required to eliminate this effect. High-precision terrain data is used to correct the terrain phase error in SAR images. The orbital data of Sentinel SAR is obtained to ensure that the SAR images acquired at different time points can be accurately aligned in space.

[0042] The GMTSAR open source software installed on the Ubuntu platform is used to perform standard differential interferometry processing: multiple acquired SAR images are georeferenced to ensure that all acquired SAR images are accurately aligned in space. Of course, before the acquired SAR images are georeferenced, conventional splicing and cropping operations are performed on the SAR images. The spatiotemporal baseline between each SAR image is calculated based on the SAR orbit data, and a reasonable spatiotemporal baseline threshold is set based on the spatiotemporal baseline. In this embodiment, the spatiotemporal baseline threshold is set to 50 days and 150 meters. Based on the spatiotemporal baseline threshold, the SBAS technology is used to select the interference combination and generate a time series interference pair.

[0043] The Goldstein filter is used to perform spatiotemporal filtering on the generated time series interference pairs, and the average coherence of the water area is collected. The average coherence of the water area is set as the threshold for interferometric phase unwrapping. The Snaphu software is used to perform phase unwrapping on the interference phase to generate the unwrapped time series interference map. The Goldstein filter is used to filter the interference map, which can reduce the impact of noise and improve image quality. The Snaphu software performs phase unwrapping on the interference pairs. During the phase unwrapping process, a coherence threshold of 0.05 is used to mask the water area in the interference map, which reduces unnecessary calculations and speeds up the unwrapping speed. The 30m resolution SRTM V3 DEM is used to remove the terrain phase in the interference map, remove the phase changes caused by the terrain from the interference map, and ensure that the final phase information only reflects the surface deformation rather than the terrain undulation.

[0044] In short, precise georeferencing, strict interferogram selection, and high-quality interferogram generation are used to ensure the accuracy and reliability of the final results.

[0045] Intermittent stacking time series InSAR processing: Check the unwrapping error and atmospheric delay of a certain pixel in the unwrapped time series interferogram, identify and remove those interferograms with missing pixels (that is, there are obvious errors or the atmospheric effect is large) to ensure the data quality of subsequent analysis. Then perform intermittent stacking time series InSAR processing on the pixel in the remaining interferograms, and calculate the annual average deformation rate of the pixel; repeat this step to calculate the annual average deformation rate of each pixel in the time series interferogram, and finally generate the annual average deformation rate result of the target area.

[0046] Specifically, in this embodiment, the principle of intermittent stacking time series InSAR processing can be found in Figure 2For the M interferograms in the time series, some interferograms (3, ..., M-1) lack effective values ​​at the pixel point (x, y). These interferograms that lack effective values ​​at the pixel point (x, y) are interferograms with obvious errors or large atmospheric effects. Therefore, only the interferograms (1, 2, ..., M-2, M) with phase values ​​at the pixel point (x, y) are considered to form a new interference pair combination m. For the pixel point (x, y), the annual average deformation rate V is calculated using the following formula LOS (x, y):

[0047]

[0048] Wherein, (x, y) represents the pixel point with row and column coordinates x and y in the geographic coordinate system, ΔΦ(x, y) represents the phase difference vector; m(x, y) is the pixel point (x, y) in the mth interference pattern in the interference pair combination; is the phase value of the pixel point (x, y) in the mth interference pattern in the interference pair combination; m(x ref ,y ref ) is the reference point (x) in the mth interference pattern of the interference pair combination. ref ,y ref ); is the reference point (x ref ,y ref ) is the phase value in the mth interference pattern in the interference pair combination; ΔT(x, y) represents the time span vector of the pixel point (x, y); ΔT m(x，y) is the time span of the pixel point (x, y) in the mth interference pattern in the interference pair combination; ∑ represents a diagonal matrix; V LOS (x, y) is the annual average deformation rate corresponding to the pixel point (x, y); <ΔΦ(x, y), ΔT(x, y)> is the vector dot product.

[0049] In this embodiment, the number of interference pair combinations of pixel points (x, y) and the coefficient matrix rank are also calculated:

[0050] Noi(x,y)=m(x,y) (3)

[0051] Roi(x,y)=rank(B(x,y)) (4)

[0052] Wherein, Noi(x, y) is the number of interference pair combinations of pixel point (x, y), and Roi(x, y) is the coefficient matrix rank of pixel point (x, y).

[0053] Repeat the calculation steps to calculate the annual average deformation rate of each pixel in the time series interference map, and finally generate the annual average deformation rate result of the target area.

[0054] The method of the present invention selects a corresponding interference pattern group for each pixel point, which can avoid the missing results caused by the intermittent coherence characteristics of the pixel points. By flexibly handling missing values, the applicability to low coherence areas is improved. Through the unwrapping error and atmospheric delay check, as well as the quality assessment of Noi and Roi, it is ensured that the interference pattern involved in the analysis has high quality and reliability. Using improved algorithms and precise mathematical models, the annual average deformation rate can be calculated more accurately to provide reliable results. This method is suitable for landslide hazard surveys in wide-area low-coherence areas, and can perform efficient and accurate landslide monitoring in a larger range.

[0055] In summary, the present invention solves the limitations of traditional methods in low-coherence environments by processing interferograms using intermittent stacking of time-series InSAR, and provides an efficient and accurate landslide monitoring method. The reliability and application value of the final deformation results are ensured through strict quality control and advanced calculation methods.

[0056] Landslide hazard generation and cataloging: Based on the annual average deformation rate results of the target area, mathematical statistics methods are used to extract the deformation area, and spatial clustering methods are used to identify and generate suspicious landslide polygons. Visual interpretation and on-site investigation are performed on the generated suspicious landslide polygons to obtain the true area of ​​the landslide polygon and catalog the landslide.

[0057] Since the error of the unwrapped phase image is mainly due to the random atmospheric turbulence error and the systematic error introduced by post-processing, it has random characteristics. According to the central limit theorem, the distribution of a large number of random variables is close to the normal distribution. For an area without ground deformation, the surface deformation results obtained by InSAR should obey the normal distribution N(0, σ), with a mathematical expectation of zero and a variance of σ 2 . However, due to the existence of noise reduction measures and random noise phase, the actual distribution obtained is close to the normal distribution N(μ, σ). For the normal distribution, the random variables distributed in the range of μ±3σ account for 99.74%. Therefore, in the embodiment of the method of the present invention, the obtained annual average deformation rate results are used to calculate the mean μ and standard deviation σ of the annual average deformation rate results. According to the principle that the error part of the unwrapped phase diagram obeys the normal distribution, the results outside the range of μ±3σ are divided into deformation areas, and the results within the range of μ±3σ are divided into non-deformation areas.

[0058] Landslides are a type of clustered point group. Therefore, the present invention uses a spatial clustering analysis method to identify suspicious landslide hazard areas on the original deformation results and generate suspicious landslide polygons. Specifically, in this embodiment, please refer to Figure 3 , the first step of spatial clustering analysis is performed based on the hot spot analysis algorithm of Getis-Ord Gi* statistics to obtain the characteristic points with high clustering. The calculation formula is as follows:

[0059]

[0060] Among them, n represents the number of all feature points, n ij represents the number of all feature points within the search distance d of the i-th feature point (i = 1, 2, ..., n; j = 1, 2, ..., n ij ; j≠i); x represents the value of the feature point; and S* represent the mean and standard deviation of all feature point values; see Figure 4 Then the cluster point calculation method is used to eliminate discrete points and extract the cluster surface. The boundaries of the cluster surface are connected to obtain the suspected landslide polygon.

[0061] Hot spot analysis is to determine whether a certain spatial element is a statistically significant hot spot by analyzing its neighboring environment. It is used to detect whether there is a clustering phenomenon in the spatial distribution of geographical phenomena. The principle of hot spot analysis is to compare the local sum of the element and its neighboring elements with the sum of all elements. When the local sum is very different from the expected local sum, a statistically significant z score can be generated, thereby identifying the clustering location of high-value or low-value elements in the spatial point group. The spatial clustering analysis method can quickly screen some deformation points with clustering and remove the interference of discrete points.

[0062] The method of the present invention can be applied to landslide deformation monitoring and landslide identification in low-coherence environments, and can solve the low coherence problem of traditional time-series InSAR technology caused by surface cover such as vegetation. At the same time, compared with the existing ISBAS, it has higher solution efficiency and can be used for landslide disaster surveys in wide-area low-coherence areas. The surface cover types in the applicable areas include but are not limited to vegetation, ice and snow types.

[0063] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field to which the present invention belongs, several simple deductions and substitutions can be made without departing from the concept of the present invention, which should be regarded as belonging to the protection scope of the present invention.

Claims

1. A surface deformation monitoring method based on intermittent stacking time series InSAR, characterized in that: The following steps are included: Data acquisition: Acquire multiple SAR images of the target area in ascending and descending orbits over a period of time, obtain high-precision terrain data of the target area, and obtain SAR orbit data; Standard differential interferometry processing: geo-register the acquired SAR images, calculate the spatiotemporal baseline between each SAR image based on the SAR orbit data, set the spatiotemporal baseline threshold based on the spatiotemporal baseline, select the interferometric combination based on the spatiotemporal baseline threshold, and generate a time series interferometric pair; perform spatiotemporal filtering and phase unwrapping on the generated time series interferometric pair to generate an unwrapped time series interferogram; Intermittent stacking time series InSAR processing: Perform unwrapping error check and atmospheric delay check on a certain pixel point in the unwrapped time series interferogram, remove the interferogram with missing pixel point, and then perform intermittent stacking time series InSAR processing on the pixel point in the remaining interferogram to calculate the annual average deformation rate of the pixel point; Repeat this step to calculate the annual average deformation rate of each pixel in the time series interference pattern, and finally generate the annual average deformation rate result of the target area; Landslide hazard generation and cataloging: Based on the annual average deformation rate results of the target area, the deformation area is extracted and the suspected landslide polygons are generated; the generated suspected landslide polygons are visually interpreted and field investigated to obtain the true area of ​​the landslide polygon and catalog the landslide.

2. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 1 is characterized in that: The GMTSAR open source software installed on the Ubuntu platform is used to perform the standard differential interferometry processing steps, and the SBAS technology is used to select the interference combination in the standard differential interferometry processing steps.

3. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 2 is characterized in that: In the standard differential interferometry processing step, the generated time series interferometry pairs are temporally and spatially filtered using a Goldstein filter.

4. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 2 is characterized in that: In the standard differential interferometry processing step, the average coherence of the water area is collected, and the average coherence of the water area is set as the threshold of the interference phase unwrapping, and the interference phase is phase unwrapped to generate the unwrapped time series interference pair.

5. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 4 is characterized in that: In the standard differential interferometry processing step, the interferometric phase was phase unwrapped using Snaphu software.

6. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 1 is characterized in that: In the intermittent stacking time series InSAR processing step, the following calculation formula is used to calculate the annual average deformation rate V of each pixel in the interference pattern: LOS (x, y): Wherein, (x, y) represents the pixel point with row and column coordinates x and y in the geographic coordinate system, ΔΦ(x, y) represents the phase difference vector; m(x, y) is the pixel point (x, y) in the mth interference pattern in the interference pair combination; is the phase value of the pixel point (x, y) in the mth interference pattern in the interference pair combination; m(x ref ,y ref ) is the reference point (x) in the mth interference pattern of the interference pair combination. ref ,y ref ); is the reference point (x ref ,y ref ) is the phase value in the mth interference pattern in the interference pair combination; ΔT(x, y) represents the time span vector of the pixel point (x, y); ΔT m (x, y) is the time span of the pixel point (x, y) in the mth interference pattern in the interference pair combination; ∑ represents the diagonal matrix, V LOS (x, y) is the annual average deformation rate corresponding to the pixel point (x, y); <ΔΦ(x, y), ΔT(x, y)> is the vector dot product.

7. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 6 is characterized in that: In the intermittent stacking time series InSAR processing step, the number of interference pair combinations and the coefficient matrix rank of each pixel are also calculated: Noi(x,y)=m(x,y) (3) Roi(x,y)=ramk(B(x,y)) (4) Wherein, Noi(x, y) is the number of interference pair combinations of pixel point (x, y), and Roi(x, y) is the coefficient matrix rank of pixel point (x, y).

8. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 1 is characterized in that: In the step of landslide hazard generation and cataloguing, mathematical statistics are used to extract deformation areas from the annual average deformation rate results, and spatial clustering methods are used to identify and generate suspicious landslide polygons.

9. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 8 is characterized in that: In the steps of landslide hazard generation and cataloguing, the steps of using mathematical statistics to extract deformation areas from the annual average deformation rate results are as follows: calculate the mean μ and standard deviation σ of the annual average deformation rate results, and according to the principle that the error part of the unwrapped phase diagram obeys the normal distribution, divide the results outside the range of μ±3σ into deformation areas, and divide the results within the range of μ±3σ into non-deformation areas.

10. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 8 is characterized in that: In the steps of landslide hazard generation and cataloguing, the steps of using the spatial clustering method to identify and generate suspicious landslide polygons based on the annual average deformation rate results are as follows: the first step of spatial clustering analysis is performed based on the hot spot analysis algorithm of the Getis-Ord Gi* statistic to obtain highly clustered feature points. The calculation formula is as follows: Among them, n represents the number of all feature points, n ij represents the number of all feature points within the search distance d of the i-th feature point (i = 1, 2, ..., n; j = 1, 2, ..., n ij ; j≠i); x represents the value of the feature point; and S* represent the mean and standard deviation of all feature point values; the aggregation point calculation method is used to eliminate discrete points and extract aggregation surfaces. The boundaries of the aggregation surfaces are connected to obtain suspicious landslide polygons.

Images