Time sequence InSAR landslide deformation monitoring fusion method based on deformation rate difference

Abstract

The invention provides a deformation rate difference-based time sequence InSAR landslide deformation monitoring fusion method, which comprises the following steps of: acquiring and preprocessing a synthetic aperture radar image of a monitoring area to obtain processed SAR image set data, then acquiring a linear deformation rate by using Stacking-InSAR, extracting PS point high-precision deformation speed by using PS-InSAR, and acquiring a nonlinear deformation rate by using SBAS-InSAR; coordinates are generated according to the deformation rate of each technical data point in combination with the orbit and geometric information, and creep, acceleration and temporary sliding deformation areas are divided by using ArcGIS and K-Means clustering; then, SBAS-InSAR processing is carried out on the newly-added image in the regional boundary, and dynamic deformation information is obtained; and finally, carrying out weighted averaging on the deformation rates of the three monitoring technologies. The landslide deformation monitoring precision under complex terrains and environments is effectively improved, errors caused by vector monitoring point density and deformation rate problems in single interference measurement are overcome, and powerful scientific support is provided by accurately monitoring and analyzing landslide evolution, determining the range, predicting the trend and early warning.

Description

Technical Field

[0001] The invention relates to a time-series InSAR landslide deformation monitoring fusion method based on deformation rate difference, and belongs to the field of landslide deformation monitoring. Background Art

[0002] Currently, the monitoring methods for landslide deformation include total station measurement, crack meter measurement, GPS measurement and synthetic aperture radar interferometry (InSAR).

[0003] Among them, total station measurement can simultaneously measure angles (horizontal angles and vertical angles) and distances. In landslide deformation monitoring, by setting up stations at control points with known coordinates, monitoring points on the landslide body are observed to obtain information such as the horizontal angle, vertical angle, and slope distance of the monitoring points. Then, based on the trigonometric function relationship, the three-dimensional coordinates of the monitoring points are calculated. After multiple measurements and by comparing the coordinate changes of the monitoring points at different times, the displacement of the monitoring points in the horizontal and vertical directions can be obtained. This method can simultaneously measure horizontal and vertical displacements and fully understand the deformation state of the landslide body, but this method requires line of sight between the monitoring point and the control point, which may be limited in areas with complex terrain, dense vegetation, or obstructions from buildings.

[0004] Crack meter measurement is to monitor the changes in the width, depth and length of cracks on the surface of the landslide body based on the crack meter. It is usually achieved by installing sensors on both sides of the crack. When the crack changes, the sensor senses the change and converts it into an electrical signal or other measurable signal. This method has high accuracy and can accurately measure small changes in parameters such as crack width. Some high-precision crack meters can measure crack changes at the millimeter level or even smaller, providing a guarantee for timely detection of landslide hazards. However, it can only monitor the changes in the cracks themselves and cannot reflect the deformation of the entire landslide body. If you rely solely on the crack meter, you may ignore the deformation of other parts of the landslide body, resulting in an incomplete assessment of the landslide hazard. In addition, its installation location is relatively limited and needs to be installed at the crack. Moreover, if the crack position changes or new cracks appear, the crack meter may need to be reinstalled or adjusted.

[0005] GPS measurement uses GPS satellite signal receiving equipment to measure the three-dimensional coordinates of the monitoring point, and obtains the slope deformation by comparing the elevation data at different times. It has the characteristics of high precision and high efficiency, and can realize real-time dynamic monitoring of large areas and long distances. It is not limited by line of sight conditions and can obtain horizontal and vertical displacement information at the same time. However, GPS signals may be affected by terrain, buildings, etc., and the accuracy will decrease in some complex environments, and the equipment cost is relatively high.

[0006] Synthetic aperture radar differential interferometry is a space-based earth observation technology that plays an important role in the fields of surface deformation monitoring. Stacking-InSAR technology is a time-series InSAR technology that uses a weighted average method to calculate the surface deformation rate. The data processing process is relatively simple and efficient, and can obtain large-scale surface deformation characteristics in a short time. Since the deformation rate is obtained only by linear stacking estimation, Stacking-InSAR technology only extracts linear deformation information. This technology has good applicability in areas with obvious linear deformation trends, but it needs to be combined with other methods for more in-depth research in complex environments. Persistent scatterer interferometry (PS-InSAR) is a technology that uses synthetic aperture radar (SAR) data to monitor surface deformation. It monitors small surface deformations by identifying point targets with stable scattering characteristics in the time series, namely persistent scatterers. PS-InSAR is suitable for large-scale and long-term ground deformation monitoring of landslides. The radar scattering characteristics remain stable for a long time and the measurement accuracy is high. However, it is difficult to monitor surface deformation in areas with complex terrain and high vegetation coverage density. Small baseline set interferometry (SBAS-InSAR) technology reduces the influence of spatiotemporal decorrelation by selecting a set of interferometric pairs of small baselines, thereby realizing the monitoring of surface deformation. Compared with PS-InSAR technology, SBAS-InSAR technology increases the temporal sampling rate and spatial density of observation data, and is not constrained by linear deformation models. It can better restore the nonlinear deformation characteristics of the surface and obtain regional time series ground deformation. However, the processing flow of SBAS-InSAR is relatively complex, the amount of calculation is large, and the selection of baselines is required to be high. If the baseline is not selected properly, errors may be introduced, affecting the accuracy of deformation monitoring. In view of the complex scattering characteristics of landslide bodies (uneven vegetation coverage) and the difference characteristics of deformation rates (measurement accuracy of different deformation rates in the same measurement area), whether it is Stacking-InSAR, PS-InSAR or SBAR-InSAR technology, a single solution method cannot solve the monitoring accuracy error caused by the difference in density of vector monitoring points caused by the complex geological conditions of the landslide, as well as the method defects of using different solution technologies based on adaptive zoning of deformation rates.

[0007] Landslides often occur under extremely complex geological conditions. Stacking-InSAR can provide deformation rate information over a large area, which can be used to preliminarily determine the approximate scope and overall deformation trend of the landslide, and provide a macroscopic deformation background. PS-InSAR has a better monitoring effect in areas with buildings or artificial facilities and areas with stable surface scattering characteristics, and can provide high-precision local deformation information of permanent scatterers. SBAS-InSAR can monitor areas with complex terrain and certain changes in vegetation coverage, and can capture detailed information on nonlinear deformation. Therefore, the integration of Stacking-InSAR, PS-InSAR and SBAS-InSAR technologies can make up for the limitations of a single technology, improve the accuracy and reliability of landslide monitoring, and has very important engineering significance for automatically selecting the appropriate time-series differential interferometry measurement method based on the slope deformation rate. Summary of the invention

[0008] To solve the above technical defects, the present invention provides a time-series InSAR landslide deformation monitoring fusion method based on deformation rate difference, which combines Stacking-InSAR, PS-InSAR and SBAS-InSAR to monitor the survey area and improve the measurement accuracy and efficiency.

[0009] Technical solution: The present invention provides a time-series InSAR landslide deformation monitoring fusion method based on deformation rate difference, comprising the following steps:

[0010] Step 1: Acquire and preprocess the synthetic aperture radar image of the monitoring area to obtain the processed SAR image set data;

[0011] Step 2: Use the interferometric stacking method to process the SAR image set data to obtain the linear deformation rate of the monitoring area; determine the interferometric stacking method high-speed deformation area of ​​the monitoring area according to the preset high-speed deformation threshold;

[0012] Step 3: Use the permanent scatterer synthetic aperture radar interferometry method to process the SAR image data set, calculate the PS deformation rate of each data point in the monitoring area, and construct the PS deformation zone based on the PS deformation rate and geographical location coordinate information of each data point;

[0013] Step 4: Use the small baseline set synthetic aperture radar interferometry method to perform interference processing on the SAR image set data to obtain the nonlinear deformation rate of each data point in the monitoring area, and construct the nonlinear deformation area based on the nonlinear deformation rate and geographical location coordinate information of the data point;

[0014] Step 5: Determine the PS accelerated deformation zone and the nonlinear accelerated deformation zone according to the PS deformation rate of the PS deformation zone, the nonlinear deformation rate of the nonlinear deformation zone, and the deformation rate threshold;

[0015] Step 6: According to the high-speed deformation zone, PS accelerated deformation zone and nonlinear accelerated deformation zone of the interference stacking method, the creep deformation zone, accelerated deformation zone and pre-slip deformation zone of the monitoring area are obtained;

[0016] Step 7: Continuously acquire synthetic aperture radar images of the monitoring area, and update the nonlinear deformation rate of the impending slip deformation zone according to the newly acquired synthetic aperture radar images of the monitoring area;

[0017] Step 8: According to the creep deformation zone, accelerated deformation zone and pre-slip deformation zone, the linear deformation rate, PS deformation rate and nonlinear deformation rate are weighted fused to obtain the time series InSAR landslide deformation monitoring fusion results.

[0018] Furthermore, in step 1, synthetic aperture radar images of the monitoring area are acquired and preprocessed to obtain processed SAR image set data, including:

[0019] Acquire M pieces of SAR image data from the Sentinel-1A satellite within the monitoring area, where M is greater than or equal to 25, perform image registration, radiation correction, and geometric correction on the acquired SAR image data, merge the processed SAR image data into set N, and obtain the processed SAR image set data.

[0020] Furthermore, in step 2, the SAR image set data is processed using the interferometric stacking method to obtain the linear deformation rate of the monitoring area, including:

[0021] The interference stacking method is used to perform stacking interference processing on the SAR image set data to obtain the landslide differential interference map;

[0022] The landslide differential interference pattern is unwrapped to obtain an unwrapped phase, and an unwrapped phase is obtained;

[0023] According to the relationship between the unwrapped phase and the surface deformation, the linear deformation rate of the area is obtained by weighted averaging.

[0024] Furthermore, according to the relationship between the unwrapped phase and the surface deformation, the method of performing weighted average to obtain the linear deformation rate of the region includes:

[0025] The average deformation rate is calculated by the formula :

[0026]

[0027] in, is the average deformation rate; n is the number of interference patterns involved in the weighted average calculation; and are the time baseline of the i-th interference pattern and its corresponding phase value respectively.

[0028] Further, according to a preset high-speed deformation threshold, the high-speed deformation zone of the monitoring area is determined by the interference stacking method, including:

[0029] The linear deformation rate and corresponding geographical location information of the monitoring area are integrated into a complete data set as the input data for cluster analysis, and the data is standardized.

[0030] The preset high-speed deformation threshold is determined based on the geological background of the study area and historical landslide monitoring data;

[0031] Through geographic information system software, the area in the monitoring area where the linear deformation rate exceeds the preset high-speed deformation threshold is determined as the interference stacking method high-speed deformation zone, and the interference stacking method high-speed deformation zone label is added to the area.

[0032] Furthermore, the SAR image data is processed using the permanent scatterer synthetic aperture radar interferometry method to calculate the PS deformation rate of the monitoring area. The PS deformation zone is constructed based on the PS deformation rate and geographic location coordinate information of the monitoring area, including:

[0033] The SAR image data is processed by using the permanent scatterer synthetic aperture radar interferometry method to calculate the PS points and PS deformation rates in the monitoring area, and the obtained PS points are screened according to the coherence threshold C to obtain the screened PS points.

[0034] The deformation rate data of permanent scatterers at the sieved PS points and the corresponding geographical location information were organized into a complete data set. The K-Means clustering method was used for cluster analysis, and clusters with similar deformation rates and similar deformation characteristics were identified as deformation clusters.

[0035] The Delaunay triangulation algorithm and geographic location coordinates are used to construct a Delaunay triangulation network for the cluster clusters in each cluster obtained after cluster analysis, and the grid is corrected according to the deformation rate and displacement direction of PS. The corrected grid boundary is the boundary of the PS deformation zone, and the PS deformation zone is obtained.

[0036] Furthermore, the coherence threshold C is set to 0.2≤C≤0.4 when the monitoring area is densely vegetated, and is set to 0.6≤C≤0.8 when the monitoring area is bare land and artificial building area.

[0037] Furthermore, the SAR image data is interferometrically processed using a small baseline set synthetic aperture radar interferometry method to obtain the nonlinear deformation rate of the monitoring area. The nonlinear deformation zone is constructed based on the nonlinear deformation rate and geographic location coordinate information of the data points, including:

[0038] Small baseline set synthetic aperture radar interferometry (SBAS-InSAR) is used to perform interference processing on the SAR image data set to obtain the differential interferogram, and the nonlinear deformation rate and deformation time series of each data point in the monitoring area are obtained by solving the problem.

[0039] The nonlinear deformation rate and corresponding geographical location information of each data point in the monitoring area are sorted into a complete data set, and the K-Means clustering method is used for cluster analysis. The clusters with similar deformation rates and similar deformation characteristics are identified as deformation clusters.

[0040] The Delaunay triangulation algorithm and geographic location coordinates are used to construct a Delaunay triangulation network for the cluster clusters in each cluster obtained after cluster analysis, and the grid is corrected according to the nonlinear deformation rate and displacement direction. The corrected grid boundary is the nonlinear deformation zone boundary, and the nonlinear deformation zone is obtained.

[0041] Further, step 5: determining the PS accelerated deformation zone and the nonlinear accelerated deformation zone according to the PS deformation rate of the PS deformation zone, the nonlinear deformation rate of the nonlinear deformation zone, and the deformation rate threshold, including:

[0042] If the PS deformation rate of the center point of each deformation zone in the PS deformation zone is greater than the deformation rate threshold I, the zone is determined to be a PS accelerated deformation zone, and a PS accelerated deformation zone label is added to the zone;

[0043] If the nonlinear deformation rate of the center point of each deformation zone of the nonlinear deformation zone is greater than the deformation rate threshold I, the zone is determined to be a nonlinear accelerated deformation zone, and a nonlinear accelerated deformation zone label is added to the zone.

[0044] Further, step 6: according to the high-speed deformation zone, PS accelerated deformation zone and nonlinear accelerated deformation zone of the interference stacking method, the creep deformation zone, accelerated deformation zone and pre-slip deformation zone of the monitoring area are obtained, including:

[0045] A region having only one of the labels of the high-speed deformation region of the interference stacking method, the PS accelerated deformation region, or the nonlinear accelerated deformation region is defined as a creep deformation region;

[0046] An area having only any two of the labels of the high-speed deformation area of ​​the interference stacking method, the PS accelerated deformation area, or the nonlinear accelerated deformation area is defined as an accelerated deformation area;

[0047] All regions with the interference stacking method high-speed deformation zone label, PS accelerated deformation zone label and nonlinear accelerated deformation zone label are defined as pre-slip deformation zones.

[0048] Further, step 7: continuously acquiring synthetic aperture radar images of the monitoring area, and updating the nonlinear deformation rate of the impending slip deformation zone according to the newly acquired synthetic aperture radar images of the monitoring area, including:

[0049] Continuously acquire synthetic aperture radar images of the monitoring area, screen the newly acquired synthetic aperture radar images of the monitoring area according to the preset coherence condition, mark the screened SAR image dataset as M1, incorporate the screened SAR image dataset M1 into the SAR image data set, and mark the updated set as N1;

[0050] The updated dataset N1 continues to use the small baseline synthetic aperture radar interferometry technique to process the SAR image data in the pre-slip deformation zone, and calculates the slope deformation rate in the pre-slip deformation zone and the signal-to-noise ratio (SNR) of SBAS-InSAR processing. B ;

[0051] The preset coherence condition is: perform image registration processing on the newly acquired SAR image, select a registered image as the main image and the coherence coefficient of each SAR image in the screened SAR image data set M1, if m1 coherence coefficients among the m coherence coefficients obtained exceed the preset threshold Q, then the registration result meets the preset coherence condition, m1>m*0.6, Q>0.6.

[0052] Further, step 8: according to the creep deformation zone, accelerated deformation zone and pre-slip deformation zone, the linear deformation rate, PS deformation rate and nonlinear deformation rate are weightedly fused to obtain the time series InSAR landslide deformation monitoring fusion results, including:

[0053] The signal-to-noise ratio (SNR) of Stacking-InSAR is calculated by statistical analysis of the signal and noise of the landslide differential interferogram. S ;

[0054] Calculation of PS-InSAR signal-to-noise ratio (SNR) using signal power and background noise power of permanent scatterers R ;

[0055] Obtaining the signal-to-noise ratio (SNR) of SBAS-InSAR processing B ;

[0056] According to the signal-to-noise ratio (SNR) of Stacking-InSAR S , PS-InSAR signal-to-noise ratio (SNR) R and the signal-to-noise ratio (SNR) of SBAS-InSAR B The linear deformation rate, PS deformation rate and nonlinear deformation rate are weightedly fused to obtain the fusion results of time-series InSAR landslide deformation monitoring.

[0057] Furthermore, the signal-to-noise ratio (SNR) of Stacking-InSAR was calculated by statistically analyzing the signal and noise of the landslide differential interferogram. S ,include:

[0058]

[0059] in, is the variance of the signal part, is the variance of the noise part,

[0060] Calculation of PS-InSAR signal-to-noise ratio (SNR) using signal power and background noise power of permanent scatterers R ,include:

[0061]

[0062] in, Signal power of permanent scatterers is the background noise power;

[0063] The signal-to-noise ratio (SNR) of the SBAS-InSAR process B The calculation methods include:

[0064]

[0065] in, is the coherence coefficient of the interference pattern, and N is the number of independent samples involved in the average.

[0066] Furthermore, according to the signal-to-noise ratio (SNR) of Stacking-InSAR S , PS-InSAR signal-to-noise ratio (SNR) R and the signal-to-noise ratio (SNR) of SBAS-InSAR B , weighted fusion of linear deformation rate, PS deformation rate and nonlinear deformation rate is performed to obtain the fusion results of time series InSAR landslide deformation monitoring, including:

[0067] The weighted average formula in the accelerated deformation zone is:

[0068] (i=1, 2, 3, j=1, 2, 3, i≠j)

[0069]

[0070]

[0071] Among them, W1, W2 and W3 are the linear deformation rate, PS deformation rate and nonlinear deformation rate obtained by Stacking-InSAR, PS-InSAR and SBAS-InSAR respectively. , and are their weighting coefficients;

[0072] The weighted average formula in the pre-slip deformation zone is:

[0073]

[0074] Where W1, W2 and W3 are the linear deformation rate, PS deformation rate and nonlinear deformation rate obtained by Stacking-InSAR, PS-InSAR and SBAS-InSAR respectively. , and are their weight coefficients, and

[0075]

[0076]

[0077]

[0078] Furthermore, the method further comprises:

[0079] Repeat steps 7 and 8 until there is no new synthetic aperture radar image of the monitoring area.

[0080] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0081] (1) The present invention fuses the results of three monitoring methods, namely Stacking-InSAR, PS-InSAR and SBAS-InSAR, to meet the measurement accuracy and efficiency requirements of different slope deformation rate zones within the same measurement area, and automatically selects the appropriate time-series differential interferometry measurement method according to the slope deformation rate, effectively improving the measurement accuracy and efficiency.

[0082] (2) Stacking-InSAR can effectively extract tiny deformation signals by stacking multi-scene interferograms to improve the signal-to-noise ratio. PS-InSAR can accurately measure highly coherent permanent scatterers and obtain high-precision deformation. SBAS-InSAR can use interferometric image pairs of small baseline subsets to obtain continuous deformation information in a large area. The combination of the three can comprehensively and accurately obtain the deformation information of the landslide at different scales, including both the overall deformation trend and local subtle changes.

[0083] (3) The combined use of Stacking-InSAR, PS-InSAR and SBAS-InSAR can complement each other and adapt to more complex monitoring environments. Stacking-InSAR can improve data utilization by fusing multi-source data and obtain effective deformation information as much as possible in different environments. PS-InSAR can play a good role in cities or areas with stable interference conditions. SBAS-InSAR has good adaptability to large-scale mountainous areas, vegetation-covered areas, etc.

[0084] (4) Stacking-InSAR can comprehensively analyze multi-period data from a holistic perspective, thereby better understanding the mechanism of landslide deformation, improving the ability to predict the development trend of landslide deformation, and providing a more scientific basis for disaster warning and prevention. Through long-term monitoring of permanent scatterers by PS-InSAR, the long-term trend of landslide deformation can be obtained; the time series analysis method of SBAS-InSAR can capture the staged changes of deformation. Joint monitoring can combine the advantages of the three technologies in time series and spatial analysis to more accurately analyze the spatiotemporal characteristics of landslide deformation. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] Figure 1 is a flow chart of the method of the present invention;

[0086] Figure 2 It is a schematic diagram of the area division of the present invention. DETAILED DESCRIPTION

[0087] The present invention will be further described below in conjunction with the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and cannot be used to limit the protection scope of the present invention.

[0088] Embodiment 1:

[0089] like Figure 1 The method flow chart of the present invention is shown as follows: Figure 1 It can be seen that the fusion method for monitoring landslide deformation based on time series InSAR with different deformation rates provided by the present invention comprises the following steps:

[0090] Step 1: Acquire and preprocess the synthetic aperture radar image of the monitoring area to obtain the processed SAR image set data;

[0091] Step 2: Use the interferometric stacking method to process the SAR image set data to obtain the linear deformation rate of the monitoring area; determine the interferometric stacking method high-speed deformation area of ​​the monitoring area according to the preset high-speed deformation threshold;

[0092] Step 3: Use the permanent scatterer synthetic aperture radar interferometry method to process the SAR image data set, calculate the PS deformation rate of each data point in the monitoring area, and construct the PS deformation zone based on the PS deformation rate and geographical location coordinate information of each data point;

[0093] Step 4: Use the small baseline set synthetic aperture radar interferometry method to perform interference processing on the SAR image set data to obtain the nonlinear deformation rate of each data point in the monitoring area, and construct the nonlinear deformation area based on the nonlinear deformation rate and geographical location coordinate information of the data point;

[0094] Step 5: Determine the PS accelerated deformation zone and the nonlinear accelerated deformation zone according to the PS deformation rate of the PS deformation zone, the nonlinear deformation rate of the nonlinear deformation zone, and the deformation rate threshold;

[0095] Step 6: According to the high-speed deformation zone, PS accelerated deformation zone and nonlinear accelerated deformation zone of the interference stacking method, the creep deformation zone, accelerated deformation zone and pre-slip deformation zone of the monitoring area are obtained;

[0096] Step 7: Continuously acquire synthetic aperture radar images of the monitoring area, and update the nonlinear deformation rate of the impending slip deformation zone according to the newly acquired synthetic aperture radar images of the monitoring area;

[0097] Step 8: According to the creep deformation zone, accelerated deformation zone and pre-slip deformation zone, the linear deformation rate, PS deformation rate and nonlinear deformation rate are weighted fused to obtain the time series InSAR landslide deformation monitoring fusion results.

[0098] Implementation principle: This embodiment fuses the results of three monitoring methods, namely Stacking-InSAR, PS-InSAR and SBAS-InSAR, to meet the measurement accuracy and efficiency requirements of different slope deformation rate zones within the same measurement area, and automatically selects the appropriate time-series differential interferometry measurement method according to the slope deformation rate, effectively improving the measurement accuracy and efficiency.

[0099] Embodiment 2:

[0100] like Figure 1 The method flow chart of the present invention is shown as follows: Figure 1 It can be seen that the fusion method for monitoring landslide deformation based on time series InSAR with different deformation rates provided by the present invention comprises the following steps:

[0101] (1) Acquire M pieces of SAR image data, where M is greater than or equal to 25, and the SAR images in the image set are image data obtained by the Sentinel-1A satellite observing the same monitoring area;

[0102] (2) Using the stacking-InSAR method, the SAR image data in step (1) are processed by stacking interferometry, where the spatial baseline is 120 m and the temporal baseline is 60 d as the threshold;

[0103] (3) The landslide differential interferogram obtained in step (2) is unwrapped to obtain the phase change, and the regional average deformation rate is obtained by weighted average according to the unwrapped phase: , the signal-to-noise ratio (SNR) of the landslide differential interferogram is calculated by statistical analysis of the signal and noise. S ;

[0104]

[0105]

[0106] in, is the average deformation rate; n is the number of interference patterns involved in the weighted average calculation; and are the time baseline of the ith interference pattern and its corresponding phase value, is the variance of the signal part obtained by Stacking-InSAR, is the variance of the noise part;

[0107] (4) The deformation rate results and corresponding geographic location information obtained by Stacking-InSAR interpretation are integrated into a complete data set as the input data for cluster analysis, and the data are standardized to eliminate the impact of different data magnitudes.

[0108] The deformation rate is obtained by phase unwrapping the differential interferogram using Stacking-InSAR, which can obtain the phase changes at different times. Based on the relationship between the phase change and the surface deformation, the average deformation rate of the surface can be calculated.

[0109] SAR images are acquired with certain geographic coordinate information, usually based on satellite orbit parameters and earth models. After the interferogram is generated, the interferogram can be geocoded to correspond to the actual geographic coordinate system.

[0110] The threshold of deformation rate is determined based on the geological background of the study area and previous monitoring data of similar landslides. The deformation rate exceeding -10 mm per year (that is, the absolute value of the deformation exceeds 10 mm. Because the deformation of the landslide is downward, it is a negative value, and its absolute value is more than 10 mm) is defined as a rapid deformation area, and the deformation rate exceeding -15 mm per year is defined as a high-speed deformation area. This part is used as the high-speed deformation area of ​​the interference stacking method for subsequent analysis.

[0111] (5) The deformation rate data obtained by Stacking-InSAR interpretation (which is the result of "performing weighted average to obtain the average deformation rate of the region" in step 3) is classified and displayed through the geographic information system (ArcGIS) software. Different colors represent different degrees of deformation. In the deformation rate map interpreted by Stacking-InSAR, the red deformation rate generally exceeds -10 mm, and different colors will be displayed in the deformation map. The high-speed deformation area is determined by the deformation rate threshold determined in step 4; the part where the deformation rate exceeds the threshold is determined as the high-speed deformation area. In the geographic information system (ArcGIS) software, red indicates areas with high deformation rate, orange indicates areas with significant deformation, and green indicates areas with low deformation rate and relative stability. The obtained high-speed deformation area (deformation rate exceeds -15 mm per year) is marked as ; The value of n represents the number of accelerated deformation zones;

[0112] (6) The SAR image data obtained in step 1 are processed using the permanent scatterer synthetic aperture radar interferometry (PS-InSAR) method to calculate the slope deformation rate in the monitoring area, and the obtained PS points are screened according to the coherence threshold C. The signal-to-noise ratio (SNR) is calculated using the signal power and background noise power of the permanent scatterer. R ;

[0113] The coherence threshold C is 0.2≤C≤0.4 when the monitoring area is densely vegetated, and 0.6≤C≤0.8 when the monitoring area is bare land and artificial building area;

[0114] In areas with dense vegetation, radar signals may be blocked or scattered by leaves and branches of plants, which will cause signal attenuation and multipath effects (i.e., signals are reflected back to the radar from different paths); these factors will reduce the coherence of the PS points. Therefore, in this environment, the coherence threshold is set lower (usually between 0.2 and 0.4) to adapt to the influence of vegetation and still screen out effective PS points. In bare land and artificial building areas, the surface is relatively flat, or artificial facilities such as buildings and roads will produce strong and stable radar reflection signals. PS points in these areas usually have higher coherence because there are no large numbers of plants or obstacles blocking the signal, so the radar signal is more easily reflected back to the radar system, and the reflected signal is more stable;

[0115] Signal-to-Noise Ratio (SNR) R for:

[0116]

[0117] in, is the signal power of the permanent scatterer, is the background noise power;

[0118] (7) K-Means clustering is an unsupervised learning algorithm that is often used to divide a data set into K different clusters. In landslide monitoring, the K-Means clustering method is used to perform cluster analysis on the deformation rate data and corresponding geographic location information obtained by PS-InSAR and SBAB-InSAR technology. The purpose is to classify pixels with similar deformation rates and geographic distribution characteristics into one category, and to divide data points into areas with more intense deformation activity and stable areas;

[0119] First, randomly select K points from the data set as the initial cluster centers. For each data point (including geographic coordinates and deformation rate), calculate its distance to all K cluster centers, the commonly used distance metric is Euclidean distance:

[0120]

[0121] in , and are the coordinates of the cluster center and the deformation rate;

[0122] Assign each data point to the cluster center that is closest to it. This is done by selecting a cluster label for each point. , so that the Euclidean distance between the point and the cluster center is minimized;

[0123]

[0124] Once all points are assigned to their corresponding clusters, the center of each cluster is recalculated. The new cluster center is the mean of all points in the cluster. For the kth cluster, its new centroid for:

[0125]

[0126] in, is the number of points in the kth cluster, is the set of all points belonging to cluster k;

[0127] Continue to iterate the above steps until the cluster center no longer changes, or the preset maximum number of iterations is reached, and each cluster after clustering is a set of point clusters obtained according to its deformation rate distribution and geographical distribution;

[0128] The deformation area is derived from the deformation cluster. After the set of pixels in a deformation cluster is Delaunay triangulated and meshed, the deformation area is obtained. The determined deformation area contains some or all of the pixels in the deformation cluster. In the K-Means clustering process, each deformation cluster has a center. This center is obtained by calculating the average value of the geographic location coordinates (usually two-dimensional or three-dimensional spatial coordinates) of all pixels in the cluster.

[0129] (8) Delaunay triangulation algorithm and geographic location coordinates are used to construct Delaunay triangulation for the PS points in each cluster obtained after cluster analysis, and the mesh is corrected according to the deformation rate and displacement direction of the PS. The corrected mesh boundary is the deformation zone boundary, and the obtained deformation zone is marked as ;

[0130] Extract the point cluster set of step (7), first select any three non-collinear points to form a triangle, then continuously add the remaining points to the triangulated network and ensure the empty circle condition. After each new point is inserted, check whether the empty circle condition is violated. If violated, adjust the triangulated network and rebuild the relevant triangles; in this way, ensure that each vertex of each triangle corresponds to a deformation point and the edges between triangles will not pass through other points. Traverse all triangles to determine whether the deformation rate of its vertex is more than twice that of the remaining vertices or the displacement direction angle is greater than 45°. If exceeded, reselect the vertex or adjust the shape of the triangle. If not exceeded, stop traversing. Compare the deformation rate of each triangle with the deformation rate of the adjacent triangles. If the deformation rates of the adjacent triangles differ greatly, their shared edges may be used as the boundary of the deformation area. According to the edges in the triangulated network, the boundaries of similar deformation areas are merged to form a continuous boundary contour. The final deformation area boundary can be represented by a polygon, which will demarcate the landslide deformation area.

[0131] The ‌empty circle condition‌ in computer graphics refers to the Delaunay empty circle rule, an algorithm for generating high-quality triangular meshes. The Delaunay empty circle rule was proposed by Russian mathematician Victor Delanoy in 1934. Its core idea is to ensure that the generated triangle circumcircle does not contain any other points‌.

[0132] (9) The SAR image data obtained in step (1) is processed using small baseline set synthetic aperture radar interferometry (SBAS-InSAR), and the deformation rate and deformation time series obtained by the solution are processed according to the K-Means clustering method and the Delaunay triangulation algorithm (such as the method in steps (7) and (8)). A Delaunay triangulation network is constructed for the PS points in each cluster obtained after cluster analysis, and the grid is corrected according to the deformation rate and displacement direction of the PS. The corrected grid boundary is the deformation zone boundary, and the obtained deformation zone is marked as , and calculate the signal-to-noise ratio (SNR) of SBAS-InSAR B :

[0133]

[0134] in, is the coherence coefficient of the interference pattern, N is the number of independent samples involved in the average;

[0135] (10) The center point of the deformation cluster in the deformation zone obtained in step (8) and step (9) is taken as the maximum deformation point, and its deformation rate is judged whether it exceeds the deformation rate threshold I, so as to determine whether the respective deformation zone is an accelerated deformation zone, and the accelerated deformation zone is recorded as PS deformation zone and nonlinear deformation zone ; and add corresponding labels in the corresponding areas.

[0136] In rock or hard soil areas (such as mountains and rocky mountainous areas), the surface hardness is usually high, the speed of landslide activity is slow, and the deformation rate is low. The deformation rate threshold I is -5mm / a-10mm / a. For soft soil or debris flow areas (such as mudstone and alluvial soil areas), landslide activities are more frequent and faster in these areas. The deformation rate threshold I is -25mm / a-15mm / a. In urban areas or artificial construction areas, due to the influence of artificial facilities such as buildings and roads, the deformation is usually more significant. The deformation rate threshold I is -10mm / a-15mm / a. In forests or vegetation areas, the radar signal is interfered by vegetation and the coherence is low. Therefore, a lower deformation rate threshold is usually required to ensure that smaller deformations can be detected. Therefore, the deformation rate threshold I is -5mm / a-15mm / a.

[0137] (11) will only occur , or The area where the three intersect (the area has only one label) is defined as the creep deformation area, the area where the two of them intersect (the area has two labels) is defined as the accelerated deformation area, and the area where the three intersect (the area has three labels) is defined as the pre-slip deformation area; the division diagram is as follows Figure 2 shown.

[0138] (12) The SAR image data obtained after step (1) are used to calculate the deformation rate and its corresponding signal-to-noise ratio (SNR) in the pre-slip deformation zone using SBAS-InSAR. B ; (Continuously acquire SAR image data: For example, the images acquired in step 1 are up to December 2023. The impact of this step is from January 2024 to the present.)

[0139] (13) If new synthetic aperture radar (SAR) image data is obtained in the monitoring area, these new data are screened. If the preset coherence condition is met, the screened SAR image data is marked as N1 and included in the SAR image data set of step (1). The updated set is recorded as N1. The screening condition is that the newly acquired SAR image data and any one of the subsequent SAR image set data in step (11) meet the preset coherence condition:

[0140] The preset coherence condition is: perform image registration processing on the newly acquired SAR image, select a registered image as the main image and the coherence coefficient of each SAR image in the newly acquired SAR image data set M1 in step (12), if m1 of the m coherence coefficients obtained exceed the predetermined threshold Q, then the registration result satisfies the preset coherence condition, m1>m*0.6, Q>0.6; registration refers to the matching of geographic coordinates of different image graphics obtained by different imaging methods in the same area, including three aspects of processing: geometric correction, projection transformation and unified scale. ‌

[0141] (14) The updated SAR image data are used to continue to use the small baseline integrated aperture radar interferometry (SBAS-InSAR) technology to calculate the deformation rate and its corresponding signal-to-noise ratio (SNR) in the pre-slip deformation zone. B ;

[0142] (15) The slope deformation results calculated by Stacking-InSAR, PS-InSAR and SBAS-InSAR in the creep deformation zone, accelerated deformation zone and pre-slip deformation zone are weighted averaged to obtain the deformation rate and time series deformation information in the deformation zone. The time series deformation information here refers to the deformation amount that increases at the monitoring point with time, which is the time-displacement curve obtained by SBAS interpretation.

[0143] If an accurate time-displacement curve is required, the time-cumulative displacement sequences calculated by Stacking-InSAR, PS-InSAR and SBAS-InSAR are obtained respectively, and the cumulative displacements at the same time point of each method are weighted according to the weighting coefficient of the following formula to obtain a new time-cumulative displacement sequence.

[0144] The weighted average formula of step 15 in the accelerated deformation zone is:

[0145] (i=1, 2, 3, j=1, 2, 3, i≠j)

[0146]

[0147]

[0148] Among them, W1, W2 and W3 are the deformation rates obtained by Stacking-InSAR, PS-InSAR and SBAS-InSAR respectively. , and are their weighting coefficients.

[0149] The weighted average formula in the pre-slip deformation zone is:

[0150]

[0151] Where W1, W2 and W3 are the deformation data obtained by Stacking-InSAR, PS-InSAR and SBAS-InSAR respectively. , and are their weight coefficients, and

[0152]

[0153]

[0154]

[0155] (16) Repeat steps (13) to (15) until there is no new SAR image data. In the landslide body, the spatial and temporal methods of landslide deformation are finally integrated to analyze the long-term trend and short-term variation of landslide deformation.

[0156] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application may adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program codes.

[0157] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0158] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0159] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0160] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A time series InSAR landslide deformation monitoring fusion method based on deformation rate difference, characterized in that: The following steps are included: Step 1: Acquire and preprocess the synthetic aperture radar image of the monitoring area to obtain the processed SAR image set data; Step 2: Use the interferometric stacking method to process the SAR image set data to obtain the linear deformation rate of the monitoring area; determine the interferometric stacking method high-speed deformation area of ​​the monitoring area according to the preset high-speed deformation threshold; Step 3: Use the permanent scatterer synthetic aperture radar interferometry method to process the SAR image data set, calculate the PS deformation rate of each data point in the monitoring area, and construct the PS deformation zone based on the PS deformation rate and geographical location coordinate information of each data point; Step 4: Use the small baseline set synthetic aperture radar interferometry method to perform interference processing on the SAR image set data to obtain the nonlinear deformation rate of each data point in the monitoring area, and construct the nonlinear deformation area based on the nonlinear deformation rate and geographical location coordinate information of the data point; Step 5: Determine the PS accelerated deformation zone and the nonlinear accelerated deformation zone according to the PS deformation rate of the PS deformation zone, the nonlinear deformation rate of the nonlinear deformation zone, and the deformation rate threshold; Step 6: According to the high-speed deformation zone, PS accelerated deformation zone and nonlinear accelerated deformation zone of the interference stacking method, the creep deformation zone, accelerated deformation zone and pre-slip deformation zone of the monitoring area are obtained; Step 7: Continuously acquire synthetic aperture radar images of the monitoring area, and update the nonlinear deformation rate of the impending slip deformation zone according to the newly acquired synthetic aperture radar images of the monitoring area; Step 8: According to the creep deformation zone, accelerated deformation zone and pre-slip deformation zone, the linear deformation rate, PS deformation rate and nonlinear deformation rate are weighted fused to obtain the time series InSAR landslide deformation monitoring fusion results.

2. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: In step 1, the synthetic aperture radar image of the monitoring area is acquired and preprocessed to obtain the processed SAR image set data, including: Obtain M pieces of SAR image data from satellites within the monitoring area, where M is greater than or equal to 25, perform image registration, radiation correction, and geometric correction on the obtained SAR image data, merge the processed SAR image data into set N, and obtain the processed SAR image set data.

3. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: In step 2, the interferometric stacking method is used to process the SAR image set data to obtain the linear deformation rate of the monitoring area, including: The interference stacking method is used to perform stacking interference processing on the SAR image set data to obtain the landslide differential interference map; Unwrapping the landslide differential interference pattern to obtain the unwrapped phase; According to the relationship between the unwrapped phase and the surface deformation, the linear deformation rate of the area is obtained by weighted averaging; According to the relationship between the unwrapped phase and the surface deformation, the method of obtaining the linear deformation rate of the area by weighted average includes: The average deformation rate is calculated by the formula : ； in, is the average deformation rate; n is the number of interference patterns involved in the weighted average calculation; and are the time baseline of the i-th interference pattern and its corresponding phase value respectively; In step 2, according to a preset high-speed deformation threshold, the high-speed deformation zone of the monitoring area using the interferometric stacking method is determined, including: The linear deformation rate and corresponding geographical location information of the monitoring area are integrated into a complete data set as the input data for cluster analysis, and the data are standardized; Through geographic information system software, the area in the monitoring area where the linear deformation rate exceeds the preset high-speed deformation threshold is determined as the interference stacking method high-speed deformation zone, and the interference stacking method high-speed deformation zone label is added to the area.

4. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: Step 3: Use the permanent scatterer synthetic aperture radar interferometry method to process the SAR image data set, calculate the PS deformation rate of the monitoring area, and construct the PS deformation zone according to the PS deformation rate and geographic location coordinate information of the monitoring area, including: The SAR image data is processed by using the permanent scatterer synthetic aperture radar interferometry method to calculate the PS points and PS deformation rates in the monitoring area, and the obtained PS points are screened according to the coherence threshold C to obtain the screened PS points. The deformation rate data of permanent scatterers at the sieved PS points and the corresponding geographical location information were organized into a complete data set. The K-Means clustering method was used for cluster analysis, and clusters with similar deformation rates and similar deformation characteristics were identified as deformation clusters. The Delaunay triangulation algorithm and geographic location coordinates are used to construct a Delaunay triangulation network for the cluster clusters in each cluster obtained after cluster analysis, and the grid is corrected according to the deformation rate and displacement direction of PS. The corrected grid boundary is the boundary of the PS deformation zone, and the PS deformation zone is obtained.

5. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: Step 4: Use the small baseline set synthetic aperture radar interferometry method to perform interference processing on the SAR image set data to obtain the nonlinear deformation rate of the monitoring area, and construct the nonlinear deformation area based on the nonlinear deformation rate and geographic location coordinate information of the data points, including: The SAR image data is processed by small baseline synthetic aperture radar interferometry to obtain differential interferograms, and the nonlinear deformation rate and deformation time series of each data point in the monitoring area are obtained by solving the problem. The nonlinear deformation rate and corresponding geographical location information of each data point in the monitoring area are sorted into a complete data set, and the K-Means clustering method is used for cluster analysis. The clusters with similar deformation rates and similar deformation characteristics are identified as deformation clusters. The Delaunay triangulation algorithm and geographic location coordinates are used to construct a Delaunay triangulation network for the cluster clusters in each cluster obtained after cluster analysis, and the grid is corrected according to the nonlinear deformation rate and displacement direction. The corrected grid boundary is the nonlinear deformation zone boundary, and the nonlinear deformation zone is obtained.

6. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: Step 5: According to the PS deformation rate of the PS deformation zone, the nonlinear deformation rate of the nonlinear deformation zone, and the deformation rate threshold, the PS accelerated deformation zone and the nonlinear accelerated deformation zone are determined, including: If the PS deformation rate of the center point of each deformation zone in the PS deformation zone is greater than the deformation rate threshold I, the zone is determined to be a PS accelerated deformation zone, and a PS accelerated deformation zone label is added to the zone; If the nonlinear deformation rate of the center point of each deformation zone of the nonlinear deformation zone is greater than the deformation rate threshold I, the zone is determined to be a nonlinear accelerated deformation zone, and a nonlinear accelerated deformation zone label is added to the zone.

7. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: Step 6: According to the high-speed deformation zone, PS accelerated deformation zone and nonlinear accelerated deformation zone of the interference stacking method, the creep deformation zone, accelerated deformation zone and pre-slip deformation zone of the monitoring area are obtained, including: A region having only one of the labels of the high-speed deformation region of the interference stacking method, the PS accelerated deformation region, or the nonlinear accelerated deformation region is defined as a creep deformation region; An area having only any two of the labels of the high-speed deformation area of ​​the interference stacking method, the PS accelerated deformation area, or the nonlinear accelerated deformation area is defined as an accelerated deformation area; All regions with the interference stacking method high-speed deformation zone label, PS accelerated deformation zone label and nonlinear accelerated deformation zone label are defined as pre-slip deformation zones.

8. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: Step 7: Continuously obtain synthetic aperture radar images of the monitoring area, and update the nonlinear deformation rate of the pre-slip deformation zone according to the newly obtained synthetic aperture radar images of the monitoring area, including: Continuously acquire synthetic aperture radar images of the monitoring area, screen the newly acquired synthetic aperture radar images of the monitoring area according to the preset coherence condition, mark the screened SAR image dataset as M1, incorporate the screened SAR image dataset M1 into the SAR image data set, and mark the updated set as N1; The updated dataset N1 continues to use the small baseline synthetic aperture radar interferometry technique to process the SAR image data in the pre-slip deformation zone, and calculates the slope deformation rate in the pre-slip deformation zone and the signal-to-noise ratio (SNR) of SBAS-InSAR processing. B ; The preset coherence condition is: perform image registration processing on the newly acquired SAR image, select a registered image as the main image and the coherence coefficient of each SAR image in the screened SAR image data set M1, if m1 coherence coefficients among the m coherence coefficients obtained exceed the preset threshold Q, then the registration result meets the preset coherence condition, wherein m1>m*0.6, Q>0.

6.

9. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: Step 8: According to the creep deformation zone, accelerated deformation zone and pre-slip deformation zone, the linear deformation rate, PS deformation rate and nonlinear deformation rate are weighted fused to obtain the time series InSAR landslide deformation monitoring fusion results, including: The signal-to-noise ratio (SNR) of Stacking-InSAR is calculated by statistical analysis of the signal and noise of the landslide differential interferogram. S ; Calculation of PS-InSAR signal-to-noise ratio (SNR) using signal power and background noise power of permanent scatterers R ; According to the signal-to-noise ratio (SNR) of Stacking-InSAR S , PS-InSAR signal-to-noise ratio (SNR) R and the signal-to-noise ratio (SNR) of SBAS-InSAR B , weighted fusion of linear deformation rate, PS deformation rate and nonlinear deformation rate is performed to obtain the fusion results of time series InSAR landslide deformation monitoring; The signal-to-noise ratio (SNR) of Stacking-InSAR is calculated by statistical analysis of the signal and noise of the landslide differential interferogram. S ,include: ； in, is the variance of the signal part, is the variance of the noise part, Calculation of PS-InSAR signal-to-noise ratio (SNR) using signal power and background noise power of permanent scatterers R ,include: ； in, Signal power of permanent scatterers, is the background noise power; The signal-to-noise ratio (SNR) of the SBAS-InSAR process B The calculation methods include: ； in, is the coherence coefficient of the interference pattern, N is the number of independent samples involved in the average; According to the signal-to-noise ratio (SNR) of Stacking-InSAR S , PS-InSAR signal-to-noise ratio (SNR) R and the signal-to-noise ratio (SNR) of SBAS-InSAR B , weighted fusion of linear deformation rate, PS deformation rate and nonlinear deformation rate is performed to obtain the fusion results of time series InSAR landslide deformation monitoring, including: The weighted average formula in the accelerated deformation zone is: i=1,2,3,j=1,2,3,i≠j; ； ； Among them, W1, W2 and W3 are the linear deformation rate, PS deformation rate and nonlinear deformation rate obtained by Stacking-InSAR, PS-InSAR and SBAS-InSAR respectively. , and are their weighting coefficients; The weighted average formula in the pre-slip deformation zone is: ； Where W1, W2 and W3 are the linear deformation rate, PS deformation rate and nonlinear deformation rate obtained by Stacking-InSAR, PS-InSAR and SBAS-InSAR respectively. , and are the weighting coefficients of linear deformation rate, PS deformation rate and nonlinear deformation rate, respectively, and ； ； 。 10. The time series InSAR landslide deformation monitoring fusion method based on deformation rate difference according to claim 1 is characterized in that: The method further comprises: Repeat steps 7 and 8 until there is no new synthetic aperture radar image of the monitoring area.

Images