A potential landslide identification method based on coupling of InSAR deformation and influencing factors

Abstract

The present invention discloses a potential landslide identification method based on the coupling of InSAR deformation and influencing factors. The method combines time-series InSAR technology, spatial statistical analysis, optical remote sensing images and machine learning methods to carry out early identification research of regional potential landslides, so as to improve the recognition efficiency and accuracy and realize the prediction of the spatial position of the landslide. The method fully considers the respective advantages of PS-InSAR and SBAS-InSAR, and extracts the surface deformation of the unit area by combining PS-InSAR and SBAS-InSAR, thereby overcoming the problem that the density of deformation points generated by SBAS-InSAR is too low and the deformation accuracy is too low in mountainous areas with more vegetation. The method comprehensively considers the formation mechanism of the landslide in the region, and on the basis of obtaining the weight of the landslide influencing factor by adopting the support vector machine model method, the landslide influencing factor is introduced into the InSAR deformation by adopting the full probability formula, and the deformation points where the influencing factor has a greater impact on the landslide are screened out, thereby improving the accuracy of the spatial clustering of the deformation points and the outlier analysis.

Description

Technical Field

[0001] The present invention relates to the field of landslide identification, and in particular to a potential landslide identification method based on the coupling of InSAR deformation and influencing factors. Background Art

[0002] Landslide refers to the geological process and natural phenomenon in which the materials that form the slope (such as rocks, sediments, artificial fills, or a combination of these materials) move downward and outward under the action of gravity and other influencing factors. Landslides not only cause damage to human life, property and buildings, seriously affecting the development of social economy, but also cause secondary disasters such as barrier lakes. Landslides are caused by a variety of factors, and there is a great deal of uncertainty about the time and place of their occurrence. Therefore, it is very important to accurately identify the location where potential landslides may occur, and how to identify potential landslides is also the key to answering the scientific question of "where are the landslides". Therefore, how to obtain the deformation of the landslide with high precision and study the monitoring technology that improves the accuracy of obtaining deformation to carry out potential landslide identification are now urgent problems in the field of landslide disasters.

[0003] Deformation monitoring technology can help people understand the deformation state of landslides and grasp the deformation stage of landslides. For landslide deformation monitoring technology, the relevant researchers initially used ground equipment (total station, crack meter, GNSS, etc.) for deformation monitoring. However, this method can only measure at a single point and cannot monitor landslide deformation over a large area. In addition, the workload is large. The environmental conditions in some landslide areas are poor and monitoring personnel are often unable to reach them. In recent years, compared with traditional landslide deformation monitoring technology, InSAR has the characteristics of wide monitoring range, high spatial resolution, good measurement accuracy, and all-day and all-weather monitoring. It has been widely used in the field of potential landslide identification, landslide deformation monitoring and early warning over a large area. So far, many researchers have used InSAR technology to obtain landslide deformation and combined with environmental factors to identify potential landslides in the study area. Limited by the principle of InSAR itself, when the deformation of the observed object is too large or the vegetation coverage rate is high, the accuracy of the obtained deformation data will be affected, resulting in the inability to identify potential landslides. Moreover, landslides are affected by both internal and external environmental factors. Minor deformations on the ground surface are not entirely caused by landslides. Therefore, it is not possible to judge potential landslides solely by relying on InSAR deformation data.

[0004] Although the research idea of ​​using multi-source data, multi-disciplinary cross-integration, and cross-border cooperation to solve the problem of potential landslide identification has received continuous attention in recent years, however, considering that most of the newly-born landslides are not among the identified potential landslide hazards, the comprehensive identification of potential landslides still needs further research. The research on the application of InSAR technology to landslide hazards mainly focuses on the deformation monitoring and deformation characteristics analysis of the known single landslides in a specific study area, and the use of this technology for the early identification of regional landslides and the determination of the scope of potential landslides. However, most of them are based on the mathematical statistical analysis of deformation point rates, and the formation of landslides is less considered to be the comprehensive influence of complex disaster-pregnant factors and external disaster-causing factors. Therefore, there are still certain problems in the identification accuracy of landslides and the determination of the scope of potential landslides. Summary of the invention

[0005] In order to overcome the shortcomings and deficiencies of the prior art, the present invention provides a potential landslide identification method based on the coupling of InSAR deformation and influencing factors. The method includes:

[0006] Step S1, collecting landslide data per unit area, including geological maps, topographic maps, high-resolution remote sensing images, road data, river data and multi-temporal optical remote sensing images;

[0007] Step S2: Based on the multi-temporal optical remote sensing images, the landslide disasters per unit area are identified by remote sensing identification methods, and the landslide information is extracted by GIS software. Meanwhile, multi-source data are integrated to extract the landslide influencing factors, analyze the engineering geological environment per unit area, and summarize the development characteristics and spatiotemporal distribution laws of landslides per unit area.

[0008] Step S3, using ArcGIS software to calculate the SAR image terrain visibility per unit area, and based on ENVI software, using PS-InSAR technology to sequentially set the amplitude deviation threshold, coherence threshold, low-pass filtering in the spatial domain and high-pass filtering in the temporal domain for the ascending orbit Sentinel-1A image covering the unit area for a preset time, to reduce the atmospheric delay phase error in the interference phase, obtain the deformation information of the stable scattering point and screen out the relatively stable PS point;

[0009] Step S4, converting the geographic coordinate system to the SAR slant range coordinate system, using the SBAS-InSAR technology as the GCP point in the SBAS processing flow to obtain the surface deformation data of the unit area in the line of sight and slope direction, and analyzing the accuracy of the deformation data through the correlation of deformation points;

[0010] Step S5, using GIS software to extract influencing indicators of landslide development, selecting evaluation factors and performing factor correlation analysis, using support vector machine to evaluate the landslide susceptibility per unit area, and obtaining the contribution value of each influencing factor to the landslide;

[0011] Step S6, using GIS software spatial analysis and support vector machine, using the full probability formula to screen the deformation points, forming high-precision, high-correlation deformation points;

[0012] Step S7, performing hot spot analysis and kernel density analysis on the selected deformation points, extracting potential landslide areas, and delineating potential landslide points using optical image features and unstable deformation thresholds;

[0013] Step S8: Verify the identification results of potential landslides in high mountain valleys by using the susceptibility area division, and perform deformation analysis and cause analysis on typical landslides in the unit area.

[0014] Furthermore, the remote sensing identification method uses the identified landslide points as samples for training machine learning, establishes a support vector machine model to classify the landslide susceptibility level per unit area, provides verification and identification means for potential landslide identification, uses the full probability formula to screen potential landslide deformation points, and uses spatial statistical analysis to extract potential landslide areas, and combines optical impact characteristics with unstable deformation thresholds to determine potential landslides.

[0015] Furthermore, the PS-InSAR technology includes:

[0016] Connection map generation: 29 connection pairs were generated from 30-scene ascending Sentinel-1A data.

[0017] Interference processing: accurately align the 29 slave images to the super master image, resample the data, and perform interference processing on the generated 29 connection pairs. Use high-precision DEM data to remove the terrain phase in the interference phase, highlight the deformation phase information, and also generate auxiliary data such as amplitude deviation index map and average intensity map;

[0018] The deformation inversion is divided into two stages. In the first stage, the residual elevation and deformation rate are estimated by using the identified permanent scatterers and linear models, and the interference pattern is re-flattened. In the second stage, the space-time filtering technology is used to remove the atmospheric phase interference, obtain the deformation information, and generate the coherence coefficient, which is used to reflect the quality of the interference phase.

[0019] Geographic coding, by setting the thresholds of amplitude deviation and coherence coefficient, obtains the deformation information of areas with high coherence and high amplitude deviation, and converts the extracted deformation information from the SAR coordinate system to the geographic coordinate system.

[0020] Furthermore, the coherence coefficient is calculated as follows:

[0021]

[0022] Among them, γ represents the coherence coefficient, M and S represent the master and slave images that make up the interference pair respectively; * represents conjugate multiplication; m and n represent the length and width of the window respectively; (i, j) represents the position of the pixel. The larger the coherence coefficient, the better the coherence, the higher the image signal-to-noise ratio, and the less affected by noise.

[0023] Furthermore, the SBAS-InSAR technology includes:

[0024] Connection diagram generation, using SBAS-InSAR technology to perform time series deformation analysis on 30 scenes of Sentinel-1A data that have been cropped and roughly registered, and invert surface deformation information;

[0025] Differential interferometry processing: When processing the time-series deformation monitoring task, all slave images are registered to the super master image, and differential interferometry processing is performed on several connection pairs generated in each set. The multi-view processing method is used to suppress the influence of speckle noise on the results. At the same time, Goldstein filtering is used to improve the clarity of the interference fringes and the phase unwrapping accuracy, and the DelaunayMCF method is used for phase unwrapping.

[0026] Track refinement and re-leveling: use the phase unwrapping result map to locate the GCP points, combine the coherence coefficient map or interference map as a reference for the selection of control points, select the GCP points, use the stable PS points obtained by combining the three thresholds of amplitude deviation, coherence coefficient and deformation rate as the GCP points for track refinement and re-leveling, and use the external DEM data to obtain the elevation information of the GCP points to eliminate the residual constant phase;

[0027] The surface deformation inversion is divided into two stages. In the first stage, the SVD method is used to combine the re-flattened interference pattern with the paired phase-elevation scaling factors to make a preliminary estimate of key information such as residual terrain and deformation. In the second stage, the obtained GCP points are used to eliminate the interference of constant phase or phase slope. At the same time, the atmospheric phase component is removed through filtering technology to improve the deformation inversion accuracy.

[0028] Geocoding: geocoding the acquired deformation information. After two inversions, the deformation information is in the SAR coordinate system. The deformation inversion result is converted from the SAR coordinate system to the geographic coordinate system to obtain the encoded deformation rate and vector deformation point.

[0029] Furthermore, the geocoding sets the thresholds of rate accuracy and elevation accuracy. The elevation accuracy is related to the coherence, baseline length and wavelength of the interference pair, and the rate accuracy is related to the coherence and wavelength of the interference pair. The calculation formulas are:

[0030]

[0031]

[0032] Among them, P H Indicates elevation accuracy, P V represents rate accuracy, γ represents coherence coefficient, λ represents wavelength, and R represents baseline length. The higher the value of elevation accuracy or rate accuracy, the lower the measurement accuracy, which will also affect the coverage of the monitoring value.

[0033] Furthermore, the total probability formula regards the occurrence of landslide as a total probability event, and calculates the probability of landslide occurrence in each influencing factor classification based on the statistical analysis results of historical landslides in different influencing factor classifications. The cumulative probability P(B) of the occurrence of landslide events in a unit area for each influencing factor is calculated based on the total probability formula, and the expression is:

[0034]

[0035] Among them, P(B) is the cumulative probability contribution of each influencing factor to the deformation point as a potential landslide point; P(A i ) is the probability contribution value of each independent influencing factor to landslide, P(B|A i ) represents the probability value of each influencing factor category under the occurrence of landslide events, i represents the number of iterations, and n represents the total number of accumulations.

[0036] Furthermore, the hotspot analysis calculates the Z score and P value of each element in the vector data set, identifies the statistically significant high-value and low-value points in the deformation points, and determines the locations where the deformation points are clustered in space. The expression is:

[0037]

[0038] in, The value x indicates the Z score. A higher Z score indicates a higher degree of aggregation of high values, and a lower Z score indicates a higher degree of aggregation of low values. j represents the attribute value of deformation point j, that is, the deformation rate value; w ij is the spatial weight value between deformation point i and deformation point j, n is the total number of deformation points, and j represents the number of iterations.

[0039] Furthermore, the kernel density analysis calculates the probability density of an element in its surrounding neighborhood according to the search radius, and the expression is:

[0040]

[0041] in, Represents the estimated density value of the (x, y) position point, based on the 4th-order kernel function, d i Indicates the distance between a point (x, y) and a point (x i ,y i ), h represents the search distance; n is the number of deformation points, i represents the number of iterations, p i Represents the search coefficient.

[0042] Beneficial effects:

[0043] The present invention proposes a potential landslide identification method based on the coupling of InSAR deformation and influencing factors. The method fully considers the respective advantages of PS-InSAR and SBAS-InSAR, and extracts the surface deformation of the unit area by combining PS-InSAR and SBAS-InSAR, thereby overcoming the problem that the density of deformation points generated by SBAS-InSAR is too low and the deformation accuracy is too low in mountainous areas with more vegetation. The method comprehensively considers the formation mechanism of the landslide in the region, and on the basis of obtaining the weight of the landslide influencing factor by adopting the support vector machine model method, the landslide influencing factor is introduced into the InSAR deformation by adopting the full probability formula, and the deformation points with greater influence of the influencing factor on the landslide are screened out, thereby improving the accuracy of spatial clustering of deformation points and outlier analysis. The method can accurately identify potential landslides and avoid the occurrence of dangerous accidents. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 is a general step flow chart of the present invention;

[0045] Figure 2 is a flow chart of the SBAS processing of the present invention;

[0046] Figure 3 The figure is a flow chart of landslide susceptibility assessment according to the present invention. DETAILED DESCRIPTION

[0047] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application may be combined with each other. The present application is further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0048] like Figure 1 As shown in Figure 1, a potential landslide identification method based on the coupling of InSAR deformation and influencing factors is proposed.

[0049] Step S1, collecting regional landslide data per unit area, including geological maps, topographic maps, high-resolution remote sensing images, road data, river data and multi-temporal optical remote sensing images.

[0050] In this embodiment, the unit area is located in Kangding County, Ganzi Prefecture, along the right bank of the Dadu River tributary. The geographical location is 101.78°-102.82° east longitude and 29.43°-30.31° north latitude. The terrain is generally high in the north and low in the south. The Yala River runs through the north and south, and the Zheduo River runs across the east and west. In addition, there are the Yulin River and Qijiagou. The overall altitude of the unit area is relatively high, spanning a wide range from 1025m to 6537m. In this vast area, most places are occupied by towering mountains, with an altitude of more than 3600m. In this vast mountain, the river is like a silver ribbon, passing through it. Along the direction of the river, the terrain is relatively low, and the altitude is mostly around 1200m. Due to long-term river scouring, deep and steep "V"-shaped river valleys have formed on both sides of the river, making the entire terrain present a unique form with low in the middle and high on both sides, forming a typical alpine canyon landform. The lithology of the strata within a unit area is complex, with strata from the Pre-Sinian to the Quaternary distributed.

[0051] The unit area is located in the subhumid climate zone of the Qinghai-Tibet Plateau and has a plateau humid climate. The climate here is dry, sunny, with a huge temperature difference between day and night, no summer all year round, and a long ice and snow period. This climatic condition gives this area a unique natural landscape and ecological environment. Near the county seat, the annual average temperature is 7.1℃, the average temperature in January is -2.5℃, and the monthly average maximum and minimum temperatures are 15.7℃ and -14.7℃ respectively, with large temperature fluctuations. The annual rainfall reaches 803.4mm, mainly concentrated from May to October, and rainfall during this period accounts for 60%-85% of the annual rainfall. Many heavy rains and continuous rains make the rainy season here particularly humid. The maximum daily rainfall is as high as 96.1mm, and the longest continuous rain can reach 58 days, with a rainfall of 427.5mm, providing development conditions for local geological disasters.

[0052] The SAR (synthetic aperture radar) image data of this embodiment is derived from the Sentinel-1 radar image data carefully developed and launched by the European Space Agency (ESA) and the European Commission for the Copernicus Earth Observation Project. It provides a rich source of information for earth science research with its characteristics of long-term, high resolution and accuracy. Sentinel-1 satellite, as a shining pearl in the field of earth observation, has a dual-constellation co-orbit polar orbit design that gives it unique observation capabilities. The dual-constellation configuration of the Sentinel-1 satellite not only ensures that the revisit period of single-constellation orbit monitoring is 12 days, but also shortens the repetition period to an amazing 6 days through the collaborative work of the dual constellations. This high-frequency observation capability enables Sentinel-1 to capture the dynamic changes of the earth's surface at an unprecedented speed.

[0053] The satellite is equipped with a C-band synthetic aperture radar based on the Prima platform, with a wavelength accurate to 0.055466 meters and a maximum coverage width of up to 400 kilometers. Sentinel-1A operates stably in a near-polar sun-synchronous orbit with an orbital inclination of 98.18° and an orbital altitude of about 693 kilometers. These parameters ensure that Sentinel-1A can provide clear, stable and high-quality radar images.

[0054] Sentinel-1A satellite provides a variety of imaging modes, including StripMap (SM), Interferometric WideSwath (IW), ExtraWideSwath (EW) and Wave (WV) imaging modes. Each mode has its own unique application scenarios and advantages. The SM mode is suitable for the management of small islands and emergency situations; the IW mode is the main choice for land observation and can meet a variety of data needs; although the EW mode has a low spatial resolution, its wide regional coverage makes it a powerful assistant for coastal monitoring; and the WV mode shows unique advantages in small-scale, high-precision imaging. In order to ensure the accuracy of SBAS-InSAR technology processing and the richness of data sources, this study selected Sentinel-1A radar satellite data in interferometric wide-swath imaging mode (IW). These data have a high resolution of 5 meters × 20 meters, a central incidence angle of 39.17°, and use VV single polarization. The original SAR image data covers almost the entire unit area, with a time span from January 3, 2021 to December 26, 2021, and a total of 30 scenes of orbit-raising data. In order to improve processing efficiency, the present invention crops the data and only retains the key information within the unit area.

[0055] DEM (digital elevation model) is a key geographic information data, which digitally simulates the surface terrain based on limited terrain elevation data, thereby realizing the digital presentation of the terrain surface morphology. In the assessment of geological disasters such as InSAR deformation monitoring and landslides, the resolution and quality of DEM data play a vital role. This is because the accuracy of DEM data directly affects the effect of terrain phase removal, and then determines the accuracy of deformation results. Therefore, selecting a reference DEM with high quality and optimal resolution is crucial for obtaining the best deformation results. In the research of this embodiment, the present invention uses the SRTM3 data with a resolution of 90m automatically downloaded by the SARscape module. SRTM data is obtained from the Space Shuttle Radar Topography Mission (SRTM) completed by the National Aeronautics and Space Administration (NASA) and the National Mapping and Mapping Agency (NIMA) of the U.S. Department of Defense. This mission made a high-precision DEM through radar images in February 2000. In the SBAS-InSAR technology processing flow, DEM is not only the basic data for obtaining terrain factors, but also the core element in the InSAR interference processing process. The present invention uses these DEM data for image registration and the elimination of terrain phase in the interference phase to ensure the accuracy of subsequent analysis. In addition, in the analysis of influencing factors of geological disasters such as landslides, DEM can not only be used to describe slope height, but also to further calculate terrain factors such as slope, slope direction, and surface undulation, providing important reference information for disaster prevention and response.

[0056] The present embodiment makes comprehensive use of a variety of optical image data in order to analyze the landslide deformation more comprehensively and meticulously. Among them, the present invention uses the Sentinel-2A data product with a resolution of 10m, which comes from the "Copernicus" program of the European Space Agency, which aims to provide high-quality earth observation data to the world through advanced satellite technology. The Sentinel-2A satellite, with its sun-synchronous orbit and the ability to revisit the same area every 10 days, provides the present invention with rich multi-spectral imaging data, which helps the present invention to deeply understand the subtle changes in the landslide terrain. In addition, the present invention also introduces the three-dimensional image data of Google Earth, with a resolution of up to 1m, and can even provide high-precision images of about 0.5m for typical buildings. Google Earth combines satellite images with aerial photography data with its powerful data integration capabilities, presenting a three-dimensional and vivid earth surface morphology for the present invention. In order to better analyze the landslide deformation, the present invention fuses these two types of data with other optical remote sensing data, making full use of their respective advantages, so as to more accurately identify the morphological characteristics of the landslide, judge the deformation trend, and provide strong support for the prevention and response of landslide disasters. The slope height and slope gradient were obtained by processing DEM using ArcGIS software.

[0057] The stratigraphic lithology data were obtained by vectorizing the 1:200000 geological map downloaded from the National Geological Data Center. The road and river data were obtained from the National Geographic Information Resource Directory Service System.

[0058] NDVI (Normalized Difference Vegetation Index) and NDBI (Normalized Difference Building Index) were obtained by performing band calculations in ENVI software on the landsat8 OLI data (30 m resolution) downloaded from the Geospatial Data Cloud (www.gscloud.cn).

[0059] Step S2: Based on multi-temporal optical remote sensing images, identify landslide hazards per unit area through remote sensing identification methods, extract landslide information using GIS software, integrate multi-source data, extract landslide influencing factors, analyze the engineering geological environment per unit area, and summarize the development characteristics and spatio-temporal distribution laws of landslides per unit area.

[0060] Step S3, using ArcGIS software to calculate the SAR image terrain visibility per unit area, and based on envi software, taking advantage of the PS-InSAR technology, the ascending orbit Sentinel-1A image covering the unit area is successively filtered through setting amplitude deviation threshold, coherence threshold, low-pass filtering in the spatial domain and high-pass filtering in the temporal domain to reduce the atmospheric delay phase error in the interference phase, obtain more reliable deformation information of stable scattering points and screen out relatively stable PS points.

[0061] InSAR technology principle, synthetic aperture radar interferometry (InSAR) technology, originated from the development of SAR technology and was deeply influenced by Young's double-slit interferometry theory. This revolutionary technology, by using two SAR images covering the same area, realizes high-precision monitoring of small deformations of the surface. In the InSAR processing process, the radar phase contained in the corresponding pixel is subtracted to generate an interference phase map, which reveals detailed information of the surface deformation for the present invention. Each pixel in the interference phase map contains rich information. Among them, the deformation phase reflects the real deformation of the surface and is the most concerned part of the present invention. The atmospheric impact phase, orbit error phase, terrain phase and noise phase are important factors affecting the quality of the interference phase map. With its unique advantages, the synthetic aperture radar sensor can emit radar waves to the ground, and realize accurate detection of the surface morphology based on the backscattering intensity of the echo signal and the phase information between the ground and the radar. In the imaging process, a coherent radar image is generated by two imagings, and these images are further conjugate multiplied to obtain an interference phase map. The interference phase diagram contains a wealth of information, which not only reflects the terrain phase of the undulating terrain, but also includes the atmospheric influence phase caused by atmospheric disturbances, the orbit error phase caused by satellite orbit errors, and the inevitable noise phase. The principle of InSAR technology is to assume that the radar satellite observes the ground target P twice at positions A1 and A2, and the distances from the satellite to P are R1 and R2 respectively; the distance between the satellite at positions A1 and A2 is B, also known as the spatial baseline; the vertical height from the satellite's first observation position A to the ground is H; the elevation of the ground target P is h; α is the angle between the spatial baseline and the horizontal line; θ is the angle of incidence. For the ground target point P, the observation phase at the radar positions A1 and A2 can be expressed as:

[0062]

[0063]

[0064] Where λ is the radar wavelength, is the random phase between the two imagings, which is generally not considered in practical applications. The difference between the two observed phases of the ground target is:

[0065]

[0066] Right now

[0067]

[0068] Since R1>>B under the conditions of satellite interferometry, the slant range difference is approximately equal to the parallel baseline, so

[0069]

[0070] Combine the above formula to obtain λ, B, θ can be used to calculate the elevation value h, and the calculation formula is as follows:

[0071]

[0072] h=Hr×cosθ

[0073] InSAR technology is based on an idealized assumption, that is, it is assumed that the target point on the ground always remains in the same position in two radar observations. However, the reality is often not the case. As time goes by, the surface target point often shifts, and this displacement may be caused by natural factors such as crustal movement, groundwater changes, or human factors such as engineering construction. Therefore, in the interference phase, in addition to the phase change caused by terrain undulations, it also inevitably includes the phase change caused by surface deformation. In order to more accurately obtain the deformation information of the surface, the present invention needs to perform differential processing on the interference pattern to remove interference factors such as terrain phase, atmospheric phase and noise phase. This differential interferometry measurement technology (D-InSAR) provides an effective means for the present invention. D-InSAR technology includes multiple methods such as dual-track method, triple-track method and quad-track method, among which the dual-track method has been widely used because of its simple operation and accurate results. The initial real interference phase It can be expressed as:

[0074]

[0075] In the formula, Indicates the surface deformation phase; represents the terrain phase; represents the reference ellipsoid phase; represents the atmospheric phase; represents the noise phase.

[0076] By removing the phases other than the deformation phase, the change in the surface deformation in the radar line of sight direction is obtained:

[0077]

[0078] Due to the complexity of radar side-view slant range imaging and ground undulation, images acquired by radar satellites usually have geometric and radiometric distortions. The geometric distortion of radar images is caused by radar side-view imaging, and the degree and pattern of distortion depend on the spatial geometric relationship between the radar line of sight and the terrain in the radar scanning area. The geometric distortions in SAR images mainly include perspective contraction, overlap and shadow. Overlap, shadow, perspective contraction and coherent speckle are several common phenomena in radar remote sensing technology, especially when processing SAR images of complex terrain. Overlap usually occurs in areas with large terrain undulations. When the radar incident angle is less than the ground slope angle, the positions of two points on the ground on the SAR image will be reversed, forming overlap. This is because the slant range of the radar beam when it reaches the high point is actually less than the slant range when it reaches the low point, resulting in the signal at the high point being received by the radar earlier, thus appearing brighter on the image. This overlap not only affects the geometric accuracy of the image, but may also lead to misjudgment of terrain information. The shadow phenomenon is caused by the inability of the radar beam to reach the obscured area. When the radar incident angle is less than the slope angle of the back slope, the radar beam cannot illuminate the back slope, causing these areas to appear dark in the SAR image. The lack of information in the shadow area constitutes an obstacle to the complete expression of the terrain. Perspective contraction is another common radar image distortion phenomenon. In undulating terrain, when the radar incident angle is greater than the slope angle, the slope surface will be compressed into a smaller imaging area on the SAR image, resulting in brighter imaging spots. The degree of this compression is related to the slope of the slope and the radar incident angle. When the slope is perpendicular to the radar incident wave, the perspective contraction phenomenon is the strongest. In addition, coherent spots are also a common radiation distortion phenomenon in radar images. Since different ground objects may exist in the same resolution unit, their scattering signals will interfere with each other to form coherent spots.

[0079] This kind of coherent speckle not only reduces the visual effect of the image, but also may affect the subsequent terrain and geomorphology analysis and interpretation work. The principle of PS-InSAR technology, PS-InSAR technology uses interference processing of multiple SAR images. From the beginning, it used the interference pattern stacking technology to extract the time series deformation, and has developed into the more representative Persistent Scatterer InSAR (PSInSAR) method with a single image as the main image and the Small Baseline Subsets (SBAS) method with multiple images as the main image. PS-InSAR technology detects PS points in SAR images and can estimate and obtain deformation information with centimeter or even millimeter measurement accuracy. The so-called PS point is a pixel point with stable phase in the time series image, which corresponds to exposed rocks, artificial buildings, etc. in the natural environment. The PS-InSAR algorithm uses a single main image for interference processing, and suppresses the influence of temporal decoherence and atmospheric delay by extracting time series coherent targets. The SBAS algorithm with multiple SAR images as the main image is mainly aimed at the situation with a limited number of images. It uses a small time baseline and space baseline for interference processing to reduce the influence of temporal decoherence and terrain on the phase.

[0080] SBAS-InSAR technology principle, small baseline interferometry (SBAS-InSAR) technology, as a further evolution of D-InSAR technology, brings a new perspective to time series analysis. This technology was originally proposed by researchers such as Berardino and Lanari in 2002, aiming to detect surface deformation more accurately. It is particularly suitable for extracting low-resolution, large-scale, and relatively slow-deformation surface deformation information, providing a powerful tool for geological environment monitoring and disaster warning. The core of SBAS-InSAR technology lies in its short temporal and spatial baseline principle. By carefully combining and generating differential interferometric pairs, this technology can effectively reduce the complexity of phase unwrapping, and accurately estimate the deformation phase and its error using the least squares method and singular value decomposition. Subsequently, through spatiotemporal filtering technology, this technology can accurately estimate and remove atmospheric phase delays, thereby revealing the true surface deformation. Compared with D-InSAR technology, SBAS-InSAR significantly reduces errors such as temporal and spatial decorrelation and atmospheric effects, and improves the accuracy and reliability of deformation measurement. Compared with PS-InSAR technology, SBAS-InSAR generates interference pairs based on multi-view main images, and has relatively low requirements for the number of images, so it is more flexible and applicable in practical applications. This feature enables SBAS-InSAR technology to perform well in various geological environments and deformation monitoring scenarios, providing valuable data support for earth science research and disaster prevention. Assuming that at a certain time (t0, t1, ..., t N) to obtain N+1 SAR images covering the study area, select one image as the main image for registration, set a suitable threshold for the spatiotemporal baseline, and combine the images smaller than the threshold to generate M interferograms, where M satisfies the following inequality:

[0081]

[0082] For A ,t B (t A <t B ) The jth interferogram generated by the SAR images acquired at two times is removed, the influence of orbit error, flat ground effect and terrain phase is removed, and the phase of high coherence points is unwrapped. The pixel coordinate (x, r) in the differential interferogram phase can be expressed as:

[0083]

[0084] In the formula, represents the deformation phase, where λ is the wavelength, and the displacement at the starting time t0 is considered to be 0, d(t A ,x,r),d(t B , x, r) represent the time t A ,t B Cumulative deformation relative to t0 in the radar line of sight Represents terrain phase represents the phase caused by atmospheric delay, Phase caused by coherent noise. The required deformation phase is obtained by removing the phase caused by terrain phase, atmospheric delay and noise. The deformation rate at different SAR acquisition times can be obtained by performing three-dimensional spatiotemporal phase unwrapping on N interference fringe patterns.

[0085] Step S4: On the basis of ensuring high coherence and good stability, the geographic coordinate system is converted to the SAR slant range coordinate system as the GCP point in the SBAS processing flow. The SBAS-InSAR technology is used to obtain the surface deformation data of the unit area in the line of sight and slope direction, and the accuracy of the deformation data is discussed through the correlation of deformation points.

[0086] PS-InSAR acquires ground control points and time-series deformation monitoring technology, which occupies an important position in current earth science research. Among many methods, PS-InSAR and SBAS-InSAR have become the two most commonly used monitoring methods due to their high efficiency and accuracy. PS-InSAR technology, namely permanent scatterer radar interferometry technology, relies on the support of multiple images and is particularly suitable for deformation monitoring in urban areas. In urban environments with tall buildings and rich surface cover, PS-InSAR can accurately identify and track those stable scatterers, providing meticulous deformation data for the present invention. However, when faced with complex terrain such as mountainous areas, high-quality PS point acquisition becomes difficult. SBAS-InSAR technology, namely small baseline subset radar interferometry technology, is good at obtaining long-term series of surface deformation information. However, this technology requires careful selection of GCP control points during use. The selection of these control points is crucial and must be kept away from deformation areas, residual terrain areas, and phase jump areas to ensure the accuracy and reliability of the monitoring results. In order to better play the advantages of the two technologies, the present invention attempts to combine PS-InSAR and SBAS-InSAR. The stable PS points extracted by PS-InSAR technology are used as ground control points in the SBAS-InSAR monitoring process, which can not only improve the accuracy of deformation monitoring, but also provide more comprehensive and in-depth surface deformation information for the present invention.

[0087] Connection diagram generation. A total of 30 connection pairs were generated from the 30-scene Sentinel-1A orbit-raising data, and the imaging time of the super main image was May 27, 2021. The absolute value of the spatial baseline is 8.37m at the minimum and 96.88m at the maximum, which is much smaller than the critical baseline. This ensures the quality of the interferometric phase and provides a solid foundation for the subsequent determination of the coherence coefficient and amplitude deviation threshold, as well as the extraction of stable PS points.

[0088] Interference processing. At this stage, the present invention first accurately aligns the 29 slave images to the super master image, and performs data resampling to ensure the accuracy and consistency of the data. Subsequently, the present invention performs interference processing on the 29 connected pairs generated, and successfully removes the terrain phase in the interference phase using high-precision DEM data, thereby highlighting the deformation phase information. In addition, the present invention also generates auxiliary data such as amplitude deviation index maps and average intensity maps, which provide an important basis for subsequent phase unwrapping and stable PS point extraction. Amplitude deviation is an important indicator for measuring the discrete characteristics of pixel radiation intensity. It effectively reflects the stability of the surface scatterer by calculating the ratio of the average radiation intensity value of the pixel in the time series at the same position to the standard deviation.

[0089]

[0090] In the formula, σA Represents the mean amplitude of the long time series, m A Represents the standard deviation of the amplitude of the long time series.

[0091] Deformation inversion. PS deformation inversion is divided into two stages. In the first stage, the present invention uses the identified "permanent scatterers" and linear models to accurately estimate the residual elevation and deformation rate, and then re-flattens the interference pattern. Entering the second stage, the present invention uses spatiotemporal filtering technology to effectively remove atmospheric phase interference, thereby obtaining more accurate deformation information. At the same time, the present invention also generates a coherence coefficient diagram, which can intuitively reflect the quality of the interference phase and can be expressed as:

[0092]

[0093] In the formula, M and S represent the master and slave images of the interference pair, respectively; * represents conjugate multiplication; m and n represent the length and width of the window, respectively; (i, j) represents the position of the pixel. The larger the coherence coefficient, the better the coherence, the higher the image signal-to-noise ratio, and the less affected by noise.

[0094] Geocoding: By setting the thresholds of amplitude deviation and coherence coefficient, deformation information of areas with high coherence and high amplitude deviation is obtained, and the extracted deformation information is converted from the SAR coordinate system to the geographic coordinate system.

[0095] After many rigorous experimental analyses, the present invention finally sets the amplitude deviation threshold to 3.5 and the coherence coefficient threshold to 0.8. This setting effectively improves the accuracy and reliability of deformation point extraction. On this basis, the present invention successfully extracted 130,678 deformation points, and the LOS deformation rates of these points are distributed in the range of -60.70 to 100.48 mm / a. In order to obtain more stable PS deformation points, the present invention further performs a threshold analysis on the deformation rate. Referring to the research results of predecessors, the present invention regards PS points with a deformation rate in the range of -0.1-0.1 mm / a as stable deformation points, and finally selects 308 stable PS points. In the process of track refinement and re-leveling, the selection of GCP is crucial. Strictly follow the GCP selection principle to ensure that the control points are evenly distributed in the study area, and at the same time eliminate PS points that are too densely distributed in the area. After screening, 225 GCP points are finally retained. Finally, these PS points are accurately converted from the geographic coordinate system to the SAR slant range coordinates.

[0096] SBAS-InSAR technology is used to monitor the surface deformation and generate connection maps. SBAS-InSAR technology generates multi-primary image interference maps based on the principle of short time-space baselines. It overcomes the influence of time-space de-correlation and can detect high-precision surface deformation information over a wide range. It is suitable for large-scale and small-scale regional surface deformation research. This embodiment uses SBAS-InSAR technology to perform time series deformation analysis on 30-scene ascending orbit Sentinel-1A data that have been pre-processed by cropping and rough registration, and inverts surface deformation information. The main process of SBAS processing is as follows: Figure 2 As shown. Connectivity graph generation processing of the input N-scene SAR images is a key step in time-series deformation monitoring. Theoretically, the maximum number of interference connection pairs can be expressed as N(N+1) / 2, which means that when the present invention has 30 descending orbit images, up to 465 interference pairs can be obtained. However, in practical applications, not all interference pairs are suitable for deformation analysis. To avoid complete loss of correlation, the present invention needs to set an appropriate spatiotemporal baseline threshold. In this processing, the present invention sets the maximum time baseline threshold to 180 days and the space baseline threshold to 5% of the critical baseline. Such settings are intended to ensure the quality and stability of the interference phase, thereby improving the accuracy of deformation monitoring. At the same time, the present invention uses a 3D unwrapping method to optimize the phase unwrapping effect and further improve the accuracy of deformation information. After screening, a total of 110 interference pairs that meet the spatiotemporal baseline thresholds were obtained from the ascending orbit data. According to the connection graph generation results of the ascending and descending orbit data, the imaging time of the super main image in the ascending orbit is May 27, 2021, and the minimum time baseline of the two orbit data is 12 days and the maximum time baseline is 180 days.

[0097] Differential interferometry processing, when processing the task of time-series deformation monitoring, it is a crucial step to align all slave images to the super master image. This step ensures the spatial and geometric consistency of all images, and provides an accurate data basis for subsequent differential interferometry processing. Subsequently, the present invention performs differential interferometry processing on several connection pairs generated in each set. In order to suppress the influence of speckle noise on the results, the present invention performs multi-view processing. This step effectively reduces the noise level and improves the reliability of the interference results. At the same time, in order to further improve the clarity of the interference fringes and the phase unwrapping accuracy, the present invention adopts Goldstein filtering processing. This method can remove part of the noise phase, significantly improve the signal-to-noise ratio of the interference pattern, and make the deformation information more prominent. In the phase unwrapping link, the present invention adopts the Delaunay MCF method. Unlike the minimum cost flow (MCF), this method uses the Delaunay triangular grid to cover those pixels above the unwrapping coherence threshold, and only performs phase unwrapping processing on high coherence pixels. This processing method not only improves the efficiency and accuracy of phase unwrapping, but also minimizes phase jumps, thereby obtaining smoother and more reliable deformation results. With the end of the interference flow process, the present invention generates a coherence coefficient graph, an interference graph, and a disentanglement result graph.

[0098] The present invention can effectively eliminate those interference image pairs with poor coherence and poor unwrapping effect. The present invention checked the coherence coefficient diagram and unwrapping diagram of 267 image pairs one by one, and after strict screening, finally eliminated 39 pairs of image pairs with unsatisfactory coherence coefficient and unwrapping results. Finally, 228 high-quality image pair results were retained for subsequent deformation inversion processing.

[0099] Track refinement and re-leveling, this process aims to effectively remove the residual terrain phase by accurately using DEM data and ground control points, thereby greatly reducing the unwrapping error, thereby ensuring the accuracy of the surface deformation inversion. When performing this step, the selection of GCP points is particularly important. In order to ensure the representativeness and uniformity of the distribution of GCP points, and to ensure that each GCP point can accurately reflect the surface deformation, in actual operation, the present invention makes full use of the phase unwrapping result map to locate the GCP points, and combines the coherence coefficient map or interference map as a reference for the selection of control points to ensure that the selected GCP points are consistent with the surface deformation characteristics and have sufficient stability. In addition, the present invention also cleverly uses the "stable PS point" obtained by combining the three thresholds of amplitude deviation, coherence coefficient and deformation rate as the GCP point for track refinement and re-leveling. This method not only improves the selection accuracy of GCP points, but also effectively avoids the subjectivity and blindness of selecting GCP points based on a single indicator. At the same time, the present invention obtains the elevation information of the GCP points with the help of external DEM data, thereby further eliminating the residual constant phase, making the deformation inversion result more accurate and reliable.

[0100] Surface deformation inversion, SBAS deformation inversion is a complex technology, usually divided into two stages, aiming to accurately reveal the surface deformation. In the first stage of inversion, the present invention mainly uses the SVD method to combine the re-flattened interference pattern with paired phase-elevation scale factors to make a preliminary estimate of key information such as residual terrain and deformation. This process lays a solid foundation for subsequent deformation analysis. The second deformation inversion is a further optimization and improvement of the results of the first stage. In this stage, the present invention makes full use of the GCP points obtained previously, eliminates the interference of constant phase or phase slope, and ensures the accuracy of deformation information. At the same time, the atmospheric phase component is removed by filtering technology to further improve the accuracy of deformation inversion. On this basis, the present invention fits the final displacement velocity model to obtain deformation information on the time series. It is worth noting that atmospheric inhomogeneity is an important factor affecting the accuracy of deformation monitoring, and each SAR image may be affected by the superposition of APS. In order to address this problem, the present invention adopts 365-day high-pass time filtering and 1000m low-pass spatial filtering to effectively remove the residual atmospheric phase in the deformation and ensure the accuracy and reliability of the deformation inversion results.

[0101] Geocoding, geocoding the acquired deformation information. Since the deformation information acquired after two inversions is in the SAR coordinate system, in order to obtain deformation information consistent with the geographic coordinates, it is necessary to convert the deformation inversion results from the SAR coordinate system to the geographic coordinate system, and obtain the encoded deformation rate and vector deformation points and other information, so as to better perform deformation analysis. When geocoding, it is necessary to set the thresholds for rate accuracy and elevation accuracy. The elevation accuracy is related to the coherence, baseline length and wavelength of the interference pair, and the rate accuracy is related to the coherence and wavelength of the interference pair. The calculation formulas are:

[0102]

[0103]

[0104] P H Indicates elevation accuracy, P V Indicates rate accuracy. The higher the value of elevation accuracy or rate accuracy, the lower the measurement accuracy, and it will also affect the coverage of the monitoring value. With the help of histogram, the accuracy information is statistically analyzed to reasonably determine the rate accuracy and elevation accuracy of the ascending track.

[0105] Step S5, using GIS to extract influencing indicators of landslide development, comprehensively considering and selecting evaluation factors and conducting factor correlation analysis, evaluating the landslide susceptibility of unit area based on support vector machine, and obtaining the contribution value of each influencing factor to landslide.

[0106] Step S6: Based on GIS spatial analysis and support vector machine, the deformation points are screened in combination with the full probability formula to form high-precision and high-correlation deformation points.

[0107] Step S7: Perform hot spot analysis and kernel density analysis on the selected deformation points to extract potential landslide areas, and delineate potential landslide points by combining optical image features and unstable deformation thresholds.

[0108] Step S8: Verify the identification results of potential landslides in high mountain valleys by using the susceptibility area division, and perform deformation analysis and cause analysis on typical landslides in the unit area.

[0109] Support Vector Machine (SVM) is a powerful machine learning algorithm, which has also shown great potential in the field of landslide research. The core idea of ​​support vector machine is to map data that is originally difficult to distinguish in low-dimensional space to high-dimensional feature space through nonlinear mapping. In this new space, the algorithm looks for an optimal classification hyperplane, which can not only classify data efficiently, but also ensure the maximization of classification intervals, thereby improving the accuracy and stability of classification. In landslide research, the present invention faces a variety of complex influencing factors, and the interaction between these factors makes the prediction and classification of landslides extremely difficult. However, by utilizing the support vector machine algorithm, the present invention can map these influencing factors to a high-dimensional feature space, and then find an optimal classification surface to effectively distinguish landslides from non-landslide areas. This not only improves the accuracy of landslide prediction, but also provides strong technical support for landslide prevention and control work.

[0110] For nonlinearly separable data (x i ,y i )(x i (∈Rd, d is the data dimension) is the evaluation factor; y i Indicates whether a landslide occurs or not; i is the number of samples). Through nonlinear mapping φ(x), the data is mapped to a certain high-order feature space. Ω*φ(x)+b=0 is the hyperplane equation. At this time, the classification interval is equal to 2 / ‖ω‖. To maximize 2 / ‖ω‖, ‖ω‖ can be made 2 Minimum, and the classification line must satisfy the constraints:

[0111] y i (ω·x i +b)≥1-ε i ,(ε i ≥0)

[0112] Where: ε i is a slack variable. While finding the classification hyperplane, ε iThe smaller the value of is, the better. Therefore, the original problem is transformed into a problem of finding ‖ω‖ under the constraints of the above formula. 2 / 2+C(∑ i ε i ) is a quadratic programming problem of the minimum value, where C is the penalty factor. The judgment function is obtained by solving:

[0113]

[0114] Where: is the kernel function. Linear kernel function, polynomial kernel function, radial basis kernel function, neural network kernel function, and hyperbolic tangent kernel function are the five commonly used kernel functions. The specific expressions are as follows:

[0115] Linear kernel function:

[0116]

[0117] Polynomial kernel function:

[0118]

[0119] Where d is the order of the polynomial kernel function.

[0120] Gaussian radial basis kernel function:

[0121]

[0122] Where d is the order of the radial basis kernel function; σ is the width of the kernel function.

[0123] Neural network kernel function:

[0124]

[0125] Where v and a are constants.

[0126] Hyperbolic tangent kernel function:

[0127] K(x i ,y j )=tanh[a(x i y j )+b]

[0128] Since the Gaussian radial basis kernel function has a good effect on both large and small samples, the Gaussian radial basis kernel function is selected to perform the processing and analysis steps in this embodiment.

[0129] like Figure 3 As shown, the evaluation process includes:

[0130] Selection of evaluation factors: Collect relevant information and conduct on-site surveys, and select ten influencing factors, including slope height, lithology, slope aspect, normalized vegetation index, river network density, road density, slope, normalized building index and terrain humidity index.

[0131] Resampling of impact factors: Grid units were selected as the basic evaluation units, and ArcGIS software was used to resample all impact factors. The grid unit size was 30m*30m.

[0132] Extraction of influencing factor attribute values: Extract attribute values ​​of training data and overall data separately.

[0133] Landslide susceptibility prediction: Construct a support vector machine model in RStudio software, use the training data to build a support vector machine landslide susceptibility evaluation model, verify its prediction accuracy with test data, select the most accurate samples that have been tested and trained to rebuild the model, input the entire data set into the reconstructed support vector machine model and calculate the landslide susceptibility index.

[0134] Susceptibility evaluation: The obtained landslide susceptibility index was imported into ArcGIS software, and the susceptibility evaluation of the study area was completed using the natural break point method based on the value.

[0135] According to the principle of the method, the susceptibility model can be divided into a deterministic model and an uncertain model including knowledge-driven and data-driven. The data-driven model can quantitatively analyze the relationship between landslides and influencing factors, and is applied to the study of landslide susceptibility evaluation in large areas because of its objectivity. With the in-depth study of susceptibility evaluation, the increase in regional area and influencing factors has brought greater challenges to the nonlinear ability of the evaluation model, and the machine learning model is widely used in landslide susceptibility analysis because of its powerful nonlinear ability. This embodiment uses the support vector machine model that is most suitable for the study area to predict landslides.

[0136] In order to obtain accurate prediction results, this embodiment selects 225 landslide points (positive samples) and 2250 random non-landslide points (negative samples) in the study area of ​​remote sensing interpretation at a ratio of 1:10 to form positive and negative samples as training samples of support vector machine. The quantitative classification results of the above influencing factor analysis are used as attribute values, and the attributes of positive and negative samples on each factor are extracted by ArcGIS software as samples for model training. 70% of the landslide and non-landslide data are randomly selected from the established sample database to form the training data of the model, and the remaining 30% are used as data for testing the model. The cross-validation method is used to determine the optimal parameter combination, and the support vector machine model is constructed through the optimal parameter combination. The overall data set is input into the model to calculate the landslide susceptibility index.

[0137] The support vector machine model can not only classify efficiently, but also deeply reveal the importance of each factor. By calculating the (Increase in node purity) of each factor, that is, the residual sum of squares, the present invention is able to quantify the impact of each variable on the heterogeneity of observations on the classification tree nodes. The larger this index is, the more significant the contribution of the variable in the classification is, thus providing a powerful tool for the present invention to compare and evaluate the importance of each variable. By comparing the "Inc Node Purity" value of each factor, it can be seen that the normalized building index and slope aspect are of significant importance, and the other factors are not much different. The order of importance is NDBI, slope aspect, elevation, surface undulation, stratum lithology, NDVI, slope, TWI, distance from the road, and distance from the river.

[0138] In order to accurately describe the geological features and landslide susceptibility of the study area, this embodiment uses a 30m grid as an evaluation unit for susceptibility evaluation. First, all grid units in the study area are converted into point features using ArcGIS software, and then the multi-value extraction to point tool is selected to extract the attribute values ​​of all grid points under each factor and import them into the trained support vector machine model to obtain the susceptibility probability value of each grid point. Finally, the grid points are converted into a grid layer again, and the susceptibility index is graded using the natural breakpoint method. This embodiment divides the entire area into five levels according to the obtained susceptibility index: low susceptibility (0-0.25), relatively low susceptibility (0.20-0.45), medium susceptibility (0.45-0.67), relatively high susceptibility (0.67-0.85) and high susceptibility (0.85-1).

[0139] In order to accurately evaluate the performance of the landslide susceptibility prediction model, the present invention needs to use evaluation indicators to comprehensively test its prediction results. In the field of landslide susceptibility evaluation, ROC curves and AUC values ​​are highly favored because of their intuitiveness and effectiveness. The ROC curve, i.e., the receiver operating characteristic curve, uses the false positive rate and the true positive rate as coordinate axes to obtain different classification results by adjusting the probability threshold. Among them, the false positive rate reflects the proportion of the model mistakenly predicting non-landslide areas as landslides, while the true positive rate reflects the ability of the model to correctly identify landslide areas. The AUC value, as the area under the ROC curve, can quantify the prediction performance of the model. The closer its value is to 1, the higher the prediction accuracy of the model. The AUC value is 0.8436, indicating that this model has good accuracy.

[0140] The total probability formula regards the occurrence of landslide as a total probability event. According to the statistical analysis results of historical landslides in different influencing factor classifications, the probability of landslide occurrence in each influencing factor classification is calculated. The cumulative probability P(B) of the occurrence of landslide events in a unit area due to each influencing factor is calculated according to the total probability formula. The expression is:

[0141]

[0142] Among them, P(B) is the cumulative probability contribution of each influencing factor to the deformation point as a potential landslide point; PA( i ) is the probability contribution value of each independent influencing factor to landslide, P(B|A i ) represents the probability value of each influencing factor category under the occurrence of landslide events, i represents the number of iterations, and n represents the total number of accumulations.

[0143] Hot spot analysis calculates the Z score and P value of each feature in the vector dataset, identifies statistically significant high and low value points in the deformation points, and determines the locations where the deformation points are clustered in space. The expression is:

[0144]

[0145] in, The value x indicates the Z score. A higher Z score indicates a higher degree of aggregation of high values, and a lower Z score indicates a higher degree of aggregation of low values. j represents the attribute value of deformation point j, that is, the deformation rate value; w ij is the spatial weight value between deformation point i and deformation point j, n is the total number of deformation points, and j represents the number of iterations.

[0146] Kernel density analysis calculates the probability density of a feature in its surrounding neighborhood based on the search radius. The expression is:

[0147]

[0148] in, Represents the estimated density value of the (x, y) position point, based on the 4th-order kernel function, d i Indicates the distance between a point (x, y) and a point (x i ,y i ), h represents the search distance; n is the number of deformation points, i represents the number of iterations, p i Represents the search coefficient.

[0149] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, the terms "set", "install", "connect", "connect", and "fix" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be indirectly connected through an intermediate medium, or it can be the internal communication of two elements. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood by specific circumstances.

[0150] Although embodiments of the present invention have been shown and described, it will be understood by those skilled in the art that various equivalent changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A potential landslide identification method based on the coupling of InSAR deformation and influencing factors, characterized by , the method comprises: Step S1, collecting landslide data per unit area, including geological maps, topographic maps, high-resolution remote sensing images, road data, river data and multi-temporal optical remote sensing images; Step S2: Based on multi-temporal optical remote sensing images, identify landslide hazards per unit area through remote sensing identification methods, extract landslide information using GIS software, integrate multi-source data, extract landslide influencing factors, analyze the engineering geological environment per unit area, and summarize the development characteristics and spatiotemporal distribution laws of landslides per unit area; Step S3, using ArcGIS software to calculate the SAR image terrain visibility per unit area, and based on envi software, using PS-InSAR technology to sequentially set the amplitude deviation threshold, coherence threshold, low-pass filtering in the spatial domain and high-pass filtering in the temporal domain to reduce the atmospheric delay phase error in the interference phase, obtain the deformation information of the stable scattering point and screen out the relatively stable PS point; Step S4, converting the geographic coordinate system to the SAR slant range coordinate system, using the SBAS-InSAR technology as the GCP point in the SBAS processing flow to obtain the surface deformation data of the unit area in the line of sight and slope direction, and analyzing the accuracy of the deformation data through the correlation of deformation points; Step S5, using GIS software to extract influencing indicators of landslide development, selecting evaluation factors and performing factor correlation analysis, using support vector machine to evaluate the landslide susceptibility per unit area, and obtaining the contribution value of each influencing factor to the landslide; Step S6, using GIS software spatial analysis and support vector machine, using the full probability formula to screen the deformation points, forming high-precision, high-correlation deformation points; Step S7, performing hot spot analysis and kernel density analysis on the selected deformation points, extracting potential landslide areas, and delineating potential landslide points using optical image features and unstable deformation thresholds; Step S8, verifying the identification results of potential landslides in high mountain valleys by using the susceptibility area division, and performing deformation analysis and cause analysis on typical landslides in the unit area; The hotspot analysis calculates the Z score and P value of each element in the vector data set, identifies the statistically significant high and low value points in the deformation points, and determines the locations where the deformation points are clustered in space. The expression is: in, The value x indicates the Z score. A higher Z score indicates a higher degree of aggregation of high values, and a lower Z score indicates a higher degree of aggregation of low values. j represents the attribute value of deformation point j, that is, the deformation rate value; w ij is the spatial weight value between deformation point i and deformation point j, n is the total number of deformation points, and j represents the number of iterations; The kernel density analysis calculates the probability density of the element in its surrounding neighborhood according to the search radius, and the expression is: in, Represents the estimated density value of the (x, y) position point, based on the 4th-order kernel function, d i Indicates the distance between a point (x, y) and a point (x i ,y i ), h represents the search distance; n is the number of deformation points, i represents the number of iterations, p i Represents the search coefficient.

2. The potential landslide identification method based on the coupling of InSAR deformation and influencing factors as claimed in claim 1 is characterized in that: The remote sensing identification method uses the identified landslide points as samples for training machine learning, establishes a support vector machine model to classify the landslide susceptibility level per unit area, provides verification and identification means for potential landslide identification, uses a full probability formula to screen potential landslide deformation points, and uses spatial statistical analysis to extract potential landslide areas, and combines optical impact characteristics with unstable deformation thresholds to determine potential landslides.

3. The potential landslide identification method based on the coupling of InSAR deformation and influencing factors as claimed in claim 1 is characterized in that: The PS-InSAR technology includes: Connection map generation: 30 connection pairs are generated from 30-scene ascending Sentinel-1A data; Interference processing: accurately align the 29 slave images to the super master image, resample the data, and perform interference processing on the generated 29 connection pairs. Use high-precision DEM data to remove the terrain phase in the interference phase, highlight the deformation phase information, and also generate auxiliary data such as amplitude deviation index map and average intensity map; The deformation inversion is divided into two stages. In the first stage, the residual elevation and deformation rate are estimated by using the identified permanent scatterers and linear models, and the interference pattern is re-flattened. In the second stage, the space-time filtering technology is used to remove the atmospheric phase interference, obtain the deformation information, and generate the coherence coefficient, which is used to reflect the quality of the interference phase. Geographic coding, by setting the thresholds of amplitude deviation and coherence coefficient, obtains the deformation information of areas with high coherence and high amplitude deviation, and converts the extracted deformation information from the SAR coordinate system to the geographic coordinate system.

4. The potential landslide identification method based on the coupling of InSAR deformation and influencing factors as claimed in claim 3 is characterized in that: The coherence coefficient is calculated as follows: Among them, γ represents the coherence coefficient, M and S represent the master and slave images that constitute the interference pair respectively; * represents conjugate multiplication; m and n represent the length and width of the window respectively; (i, j) represents the position of the pixel. The larger the coherence coefficient, the better the coherence, the higher the image signal-to-noise ratio, and the less affected by noise.

5. The potential landslide identification method based on the coupling of InSAR deformation and influencing factors as claimed in claim 1 is characterized in that: The SBAS-InSAR technology includes: Connection diagram generation, using SBAS-InSAR technology to perform time series deformation analysis on 30 scenes of Sentinel-1A data that have been cropped and roughly registered, and invert surface deformation information; Differential interferometry processing: When processing the time-series deformation monitoring task, all slave images are registered to the super master image, and differential interferometry processing is performed on several connection pairs generated in each set. The multi-view processing method is used to suppress the influence of speckle noise on the results. At the same time, Goldstein filtering is used to improve the clarity of the interference fringes and the phase unwrapping accuracy, and the DelaunayMCF method is used for phase unwrapping. Track refinement and re-leveling: use the phase unwrapping result map to locate the GCP points, combine the coherence coefficient map or interference map as a reference for the selection of control points, select the GCP points, use the stable PS points obtained by combining the three thresholds of amplitude deviation, coherence coefficient and deformation rate as the GCP points for track refinement and re-leveling, and use the external DEM data to obtain the elevation information of the GCP points to eliminate the residual constant phase; The surface deformation inversion is divided into two stages. In the first stage, the SVD method is used to combine the re-flattened interference pattern with the paired phase-elevation scaling factors to make a preliminary estimate of the residual terrain and deformation information. In the second stage, the obtained GCP points are used to eliminate the interference of constant phase or phase slope. At the same time, the atmospheric phase component is removed through filtering technology to improve the deformation inversion accuracy. Geocoding: geocoding the acquired deformation information. After two inversions, the deformation information is in the SAR coordinate system. The deformation inversion result is converted from the SAR coordinate system to the geographic coordinate system to obtain the encoded deformation rate and vector deformation point.

6. A potential landslide identification method based on coupling of InSAR deformation and influencing factors as claimed in claim 5, characterized in that: The geocoding is used to set the thresholds of rate accuracy and elevation accuracy. The elevation accuracy is related to the coherence, baseline length and wavelength of the interference pair, and the rate accuracy is related to the coherence and wavelength of the interference pair. The calculation formulas are: Among them, P H Indicates elevation accuracy, P V represents rate accuracy, γ represents coherence coefficient, λ represents wavelength, and R represents baseline length. The higher the value of elevation accuracy or rate accuracy, the lower the measurement accuracy, which will also affect the coverage of the monitoring value.

7. The potential landslide identification method based on the coupling of InSAR deformation and influencing factors as claimed in claim 1 is characterized in that: The total probability formula regards the occurrence of landslide as a total probability event. According to the statistical analysis results of historical landslides in different influencing factor classifications, the probability of landslide occurrence in each influencing factor classification is calculated. The cumulative probability P(B) of the occurrence of landslide events for each influencing factor in a unit area is calculated according to the total probability formula. The expression is: Among them, P(B) is the cumulative probability contribution of each influencing factor to the deformation point as a potential landslide point; P(A i ) is the probability contribution value of each independent influencing factor to landslide, P(B|A i ) represents the probability value of each influencing factor category under the occurrence of landslide events, i represents the number of iterations, and n represents the total number of accumulations.

Images