A surface deformation monitoring method based on intermittent stacking time series InSAR

Abstract

The present invention belongs to the field of terrain remote sensing monitoring technology, and specifically relates to a surface deformation monitoring method based on intermittent stacking time-series InSAR. The method first obtains multiple SAR images, high-precision terrain data, and orbit data of the target area over a period of time. The SAR images are then georeferenced, and an interference combination is selected. After spatiotemporal filtering and phase unwrapping, a time series interferogram is generated. The annual average deformation rate of each pixel in the interferogram is then calculated to generate the annual average deformation rate result of the target area. Finally, a suspected landslide polygon is generated based on the annual average deformation rate result, and visual interpretation and on-site investigation are performed to obtain the true area of ​​the landslide polygon. This method can avoid the loss of results caused by the intermittent coherence characteristics of pixel points, has a higher tolerance in pixel selection, can achieve higher coverage, significantly reduces the demand for computing power, and can quickly obtain deformation rate results in low-coherence environments.

Description

Technical Field

[0001] The present 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 has applications 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 is effective for monitoring surface deformation. However, due to the steep terrain and dense vegetation in river canyons, shortwave SAR satellites like Sentinel-1 have very low coherence in these areas. This is especially true for long-term landslide monitoring scenarios, where the acquired images exhibit intermittent coherence (intermittent coherence refers to the difficulty of ground scatterers maintaining high coherence over long time series). As a result, traditional time-series InSAR methods are able to obtain very few deformation data in low-coherence areas, making them ineffective for effective landslide monitoring and identification.

[0004] Stacking-InSAR (Stacked Synthetic Aperture Radar Interferometry) and SBAS-InSAR (Small Baseline Subset Synthetic Aperture Radar Interferometry) are commonly used time-series InSAR techniques in the field of geological disaster deformation monitoring. In 1998, Sandwell et al. proposed the Stacking-InSAR algorithm. By combining M unwrapped interferometric phase images and performing a weighted superposition, the algorithm reduces the influence of orbital, atmospheric, and topographic 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 interferometric 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 each interference pair by setting a spatiotemporal baseline threshold. At the same time, the influence of atmospheric phase can be weakened by using high-pass filtering in the time domain and low-pass filtering in the spatial domain.

[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 over 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 linear velocity of fault deformation. In 2013, Sowter et al. introduced the intermittent coherence method based on the SBAS algorithm, increasing the SBAS solution scale from the entire image to the pixel level. The proposed ISBAS algorithm can achieve a coverage range 4-40 times that of the SBAS algorithm. However, in the actual needs of deformation monitoring and geological disaster monitoring, the SBAS algorithm has a long computation time and high computer memory requirements, 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 of a set of unwrapped interferograms after quality screening in an intermittent manner. However, the scheme screens the entire interference image set according to certain standards and then manually 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 interferograms. 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 mainly targets a certain moment in the earthquake time, 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, and does not process 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 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. This method selects a corresponding interferogram for each pixel, avoiding missing results caused by the intermittent coherence characteristics of the pixels. In pixel selection, missing values ​​are allowed in the interferograms at some time points, and only those with phase values ​​are considered. Compared with existing ISBAS technology, this method has higher tolerance and can achieve higher coverage. By flexibly handling missing values ​​and low-coherence areas, this method significantly reduces the demand for computing power, enabling rapid acquisition of InSAR deformation rate results in low-coherence environments, and realizing rapid landslide hazard surveys in wide-area low-coherence environments. It has great application value in the field of landslide and other geological disaster monitoring.

[0014] The method includes the following steps: data acquisition: acquiring multiple SAR images of the target area in ascending and descending orbits over a period of time, acquiring high-precision terrain data of the target area, and acquiring SAR orbit data. Standard differential interferometry processing: georeferencing the acquired SAR images, calculating the spatiotemporal baseline between each SAR image based on the SAR orbit data, setting the spatiotemporal baseline threshold based on the spatiotemporal baseline, selecting an 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 an unwrapping error and atmospheric delay check on a certain pixel point in the unwrapped time series interference map, eliminating the interference map with the pixel point missing, and then performing intermittent stacking time series InSAR processing on the pixel point in the remaining interference map to calculate the annual average deformation rate of the pixel point; repeating this step to calculate 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 suspicious landslide polygons are generated; the generated suspicious landslide polygons are visually interpreted and field surveyed to obtain the true area of ​​the landslide polygons and conduct landslide cataloging.

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

[0016] In a specific embodiment, 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 embodiment, 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, during the standard differential interferometry processing step, Snaphu software is used to perform phase unwrapping on the interferometric phase.

[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 interferogram: LOS (x, y):

[0020]

[0021] Where (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 interferogram of 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 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.

[0022] In a specific embodiment, 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 embodiment, in the step of generating and cataloguing landslide hazards, the deformation area is extracted by applying mathematical statistics to 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 embodiment, in the step of generating and cataloging landslide hazards, the step of using mathematical statistics to extract the deformation area of ​​the annual average deformation rate results is: calculating 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, 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.

[0028] In a specific embodiment, in the step of generating and cataloging landslide hazards, the step of using a spatial clustering method to identify and generate suspicious landslide polygons based on the annual average deformation rate results is as follows: a first step of spatial clustering analysis is performed based on a hot spot analysis algorithm based on 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 Indicates 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 element 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 the aggregation surface. The boundaries of the aggregation surface are connected to obtain the suspicious landslide polygon.

[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 handling 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. It 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 Flowchart 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 Schematic diagram of the process from extracting to aggregating surfaces to generating suspicious 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 acquire multiple SAR images of the target area over a period of time, both ascending and descending, to obtain high-precision terrain data of the target area, as well as 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 to obtain more comprehensive information on surface deformation, images from different orbital directions can provide complementary information, helping to improve the quality of the interferogram and the accuracy of the interpretation. High-precision terrain data refers to digital elevation models (DEMs), which provide detailed terrain information within the study area, including altitude, slope, etc. Because the SAR signal propagation path is affected by terrain undulations, an accurate DEM is required to eliminate this effect. High-precision terrain data is used to correct terrain phase errors in SAR images. The orbital data of the Sentinel SAR is obtained to ensure that SAR images acquired at different time points are accurately aligned in space.

[0042] GMTSAR open-source software, installed on an Ubuntu platform, performs standard differential interferometry processing. Multiple acquired SAR images are georeferenced to ensure precise spatial alignment. Conventional SAR image stitching and cropping are also performed before georeferencing. The spatiotemporal baselines between the SAR images are calculated based on SAR orbital data. A reasonable spatiotemporal baseline threshold is set based on this spatiotemporal baseline. In this example, the spatiotemporal baseline thresholds are set to 50 days and 150 meters. Based on these spatiotemporal baseline thresholds, SBAS technology is used to select interferometric combinations and generate time series interferometric pairs.

[0043] The generated time series interferograms were spatiotemporally filtered using a Goldstein filter. The average coherence of the water region was collected and used as the threshold for interferometric phase unwrapping. Snaphu software was then used to perform phase unwrapping on the interferograms, generating an unwrapped time series interferogram. Filtering the interferograms with a Goldstein filter can reduce the impact of noise and improve image quality. Snaphu software was used to perform phase unwrapping on the interferograms. During the phase unwrapping process, a coherence threshold of 0.05 was used to mask the water regions in the interferograms, reducing unnecessary computational effort and accelerating unwrapping. A 30m-resolution SRTM V3 DEM was used to remove the terrain phase from the interferograms. Phase variations caused by topography were removed from the interferograms, ensuring that the resulting phase information reflects only surface deformation, not topographic fluctuations.

[0044] In short, precise georeferencing, strict interferometric pair 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: A pixel in the unwrapped time series interferogram is checked for unwrapping error and atmospheric delay. Interferograms with missing pixels (i.e., significant errors or significant atmospheric effects) are identified and removed to ensure data quality for subsequent analysis. Intermittent Stacking Time Series InSAR processing is then performed on the remaining interferograms for that pixel, and the annual average deformation rate for that pixel is calculated. This step is repeated to calculate the annual average deformation rate for each pixel in the time series interferogram, ultimately generating the annual average deformation rate for 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) are missing valid values ​​at the pixel point (x, y). These interferograms with missing valid 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 interferogram 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] Where (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 interferogram of 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 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.

[0049] In this embodiment, the number of interference pair combinations and the coefficient matrix rank of the pixel point (x, y) 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 this calculation 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.

[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. By checking the unwrapping error and atmospheric delay, as well as the quality assessment of Noi and Roi, it is ensured that the interference patterns involved in the analysis have high quality and reliability. Using improved algorithms and precise mathematical models, the annual average deformation rate can be calculated more accurately, providing reliable results. This method is suitable for landslide hazard surveys in wide-area low-coherence areas and can perform efficient and accurate landslide monitoring over a larger area.

[0055] In summary, this method, through intermittent stacking of time-series InSAR interferograms, overcomes the limitations of traditional methods in low-coherence environments and provides a highly efficient and accurate means of landslide monitoring. Through rigorous quality control and advanced computational methods, the reliability and application value of the final deformation results are ensured.

[0056] Landslide Hazard Generation and Cataloging: Based on the annual average deformation rate of the target area, mathematical statistics are used to extract deformation areas. Spatial clustering methods are then used to identify and generate suspected landslide polygons. Visual interpretation and field investigation of the generated suspected landslide polygons are conducted to determine the true area of ​​the landslide polygons and to catalog the landslides.

[0057] Since the error of the unwrapped phase image is mainly due to random atmospheric turbulence errors and systematic errors 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 the deformation area, and the results within the range of μ±3σ are divided into the non-deformation area.

[0058] Landslide is a type of clustered point group. Therefore, the embodiment of 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 , based on the hot spot analysis algorithm of Getis-Ord Gi* statistic, the first step of spatial cluster analysis is performed 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 Indicates 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 element point; and S* represent the mean and standard deviation of all feature point values; see Figure 4 Then, the clustering point calculation method is used to eliminate discrete points and extract the clustering surface. The boundaries of the clustering surface are connected to obtain the suspicious landslide polygon.

[0061] Hotspot analysis analyzes the neighborhood of a spatial feature to determine whether it is a statistically significant hotspot. It is used to detect clustering in the spatial distribution of geographic phenomena. Hotspot analysis works by comparing the local sum of a feature and its neighboring features with the sum of all features. When the local sum differs significantly from the expected local sum, a statistically significant z-score is generated, identifying clusters of high-value or low-value features within a spatial point cluster. Spatial cluster analysis methods can quickly screen for clustered deformation points and remove interference from 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 and ice and snow types.

[0063] The above content is a further detailed description of the present invention in conjunction with specific preferred embodiments, and the specific implementation of the present invention should not be considered to be limited to these descriptions. For those skilled in the art of the present invention, several simple deductions and substitutions can be made without departing from the concept of the present invention, and all of these should be considered to fall within the scope of protection 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 during 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: georegister 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 interferometry combination based on the spatiotemporal baseline threshold, and generate a time series interferometry pair; perform spatiotemporal filtering and phase unwrapping on the generated time series interferometry pair to generate the unwrapped time series interferogram; Intermittent stacking time series InSAR processing: Perform unwrapping error and atmospheric delay checks on a certain pixel in the unwrapped time series interferogram, remove the interferogram with the missing pixel, and then perform intermittent stacking time series InSAR processing on the pixel in the remaining interferograms to 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; Landslide hazard generation and cataloging: Based on the annual average deformation rate results of the target area, deformation areas are extracted and suspected landslide polygons are generated. Visual interpretation and on-site investigation of the generated suspected landslide polygons are performed to obtain the true area of ​​the landslide polygons and conduct landslide cataloging. 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 interferogram LOS (x, y): Among them, (x, y) represents the pixel point with row and column coordinates x and y respectively in the geographic coordinate system. represents the phase difference vector; m(x, y) is the pixel point (x, y) in the mth interference pattern of the interference pair combination; (x, y) 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 ref ,y ref ); (x ref ,y ref ) is the reference point (x ref ,y ref ) the phase value in the mth interferogram in the interference pair combination; A time span vector representing a pixel point (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); is the vector dot product.

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 set as the threshold for interferometric phase unwrapping. The interferometric phase is phase unwrapped to generate an unwrapped time series interferometric 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 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)=rank(B(x,y)) (4) Among them, 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).

7. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 1 is characterized in that: In the landslide hazard generation and cataloging steps, the deformation area is extracted by using mathematical statistics on the annual average deformation rate results, and the spatial clustering method is used to identify and generate suspicious landslide polygons.

8. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 7 is characterized in that: In the steps of landslide hazard generation and cataloging, the steps of using mathematical statistics to extract the deformation area of ​​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 map obeys the normal distribution, divide the results outside the range of μ±3σ into the deformation area, and divide the results within the range of μ±3σ into the non-deformation area.

9. The surface deformation monitoring method based on intermittent stacking time series InSAR according to claim 7 is characterized in that: In the landslide hazard generation and cataloging steps, the spatial clustering method is used to identify and generate suspicious landslide polygons based on the annual average deformation rate results. The first step is to perform spatial cluster analysis 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 Indicates 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 element 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 the aggregation surface. The boundaries of the aggregation surface are connected to obtain the suspicious landslide polygon.

Images