A landslide identification method, system, and medium

Abstract

The application discloses a landslide identification method and system and a medium, and belongs to the technical field of geological disaster monitoring and remote sensing image processing; the method comprises the following steps: acquiring a phase gradient graph and a deformation rate graph of a target region; calculating a local adaptive threshold based on the phase gradient graph, and classifying pixels into landslide pixels and background pixels by comparison; generating a geometric boundary box covering each landslide region based on the landslide pixels; extracting a corresponding deformation rate distribution graph by using each boundary box; determining a segmentation threshold based on the gray histogram features of each distribution graph, and extracting a significant deformation region according to the segmentation threshold; performing morphological processing and merging on each deformation region to obtain a result graph indicating the spatial distribution of the landslide region; and the application realizes automatic, high-precision identification and boundary positioning of large-range landslide hazards.

Description

Technical Field

[0001] This invention relates to the field of geological disaster monitoring and remote sensing image processing technology, and in particular to a landslide identification method, system and medium. Background Technology

[0002] Landslides, as a common geological hazard, are characterized by their suddenness, destructive power, and wide impact, seriously threatening people's lives and property and the stable operation of major projects. Timely and accurate identification of landslide hazards is of great significance for disaster early warning, risk assessment, and the deployment of prevention and control projects.

[0003] In recent years, Interferometric Synthetic Aperture Radar (InSAR) technology has become an important technical means for large-scale surface deformation monitoring and landslide hazard identification due to its advantages such as wide coverage, high monitoring accuracy, ability to penetrate clouds and fog, and all-weather observation. Currently, InSAR-based landslide identification methods mainly include the following categories: Manual identification and selection: combining InSAR deformation results with optical image features for manual interpretation. This method has high accuracy, but is time-consuming, labor-intensive, inefficient, and highly dependent on the professional level of technicians, making it difficult to apply to large-scale regional surveys. Automatic segmentation methods based on a single data source: for example, directly segmenting the InSAR deformation rate map using global or fixed thresholds, or using phase gradient maps (such as the results of improved phase gradient overlay methods) for binarization extraction. While these methods achieve automated processing, they have significant limitations: segmentation methods based on deformation rate maps struggle to adapt to the spatial heterogeneity of deformation rates within different landslides, easily leading to missed or false detections, and the extracted boundaries are coarse; methods relying solely on phase gradient maps often produce preliminary identification results with incomplete and discontinuous boundaries, potentially containing noise, and failing to accurately depict the complete morphology of the landslide. Machine learning and deep learning methods: While automatically identifying landslides by building models possesses some automation capabilities, they typically require a large number of high-quality training samples, the model training process is complex, computationally expensive, and prone to overfitting, limiting generalization ability and making them difficult to directly apply to large areas lacking labeled data.

[0004] Therefore, existing technologies generally suffer from the problem of balancing automation and recognition accuracy: highly efficient fully automated methods often have inaccurate boundary extraction, while methods that can obtain relatively accurate boundaries (such as manual interpretation) are inefficient. There is an urgent need for a technical solution that can efficiently, automatically, and accurately identify and locate the boundaries of large-scale landslide hazards. Summary of the Invention

[0005] The purpose of this invention is to provide a landslide identification method that enables automated, high-precision identification and boundary positioning of large-scale landslide hazards.

[0006] To address the aforementioned technical problems, embodiments of the present invention provide a landslide identification method, comprising the following steps: Acquire the phase gradient map and deformation rate map of the target region within a preset time period; Based on the features of neighboring pixels and the pixel value of each pixel in the phase gradient map, a corresponding local adaptive threshold is calculated for each pixel. By comparing the pixel value of each pixel with the local adaptive threshold of the pixel itself, the pixels are classified into landslide area pixels and background area pixels. Based on the pixels of the landslide area, generate multiple geometric bounding boxes that cover different landslide areas; The corresponding deformation rate distribution map is extracted from the deformation rate map using each geometric bounding box; Based on the grayscale histogram features of each deformation rate distribution map, the segmentation threshold of each deformation rate distribution map is determined; the corresponding deformation rate distribution map is segmented according to the segmentation threshold, and the deformation regions in the deformation rate distribution map that are higher than the segmentation threshold are extracted. Morphological processing and merging of each deformation region yielded a map showing the spatial distribution of the landslide area.

[0007] In some optional embodiments, obtaining the phase gradient map includes the following steps: A preset number of SAR images are registered to the same master image. Multiple interferometric combinations are generated based on the spatiotemporal baseline threshold. The wrapped phase representation of each interferogram is obtained through interferometric processing: ; In the formula, The winding phase in the i-th interferogram. x and y These represent the number of rows and columns of a pixel in the interferogram, respectively. For terrain phase, For deformation phase, For atmospheric delayed phase, Noise phase; From the interferograms, select M interferograms, and for each interferogram, calculate its phase gradient along at least four directions, where the phase gradient along the southeast direction with step size s is calculated as follows: ; In the formula, For the i-th interferogram at pixel The phase gradient value along the southeast (SE) direction; and For the i-th interferogram in coordinates and The winding phase value at the specified point; s is the step size used when calculating the gradient. The phase gradients of the M interferograms in the same direction are superimposed in the time dimension to obtain the stacked phase gradient in that direction: ; In the formula, For in pixels At this location, the phase gradient value after stacking along the southeast direction; M is the number of high-quality interferograms processed by phase gradient stacking. Median filtering is applied to the stacked phase gradients of all pixels to obtain the filtered phase gradient map: ; In the formula, This is the phase gradient map after median filtering along the southeast direction; (.) represents the median filter function; Repeat the above steps to obtain filtered phase gradient maps in at least four directions, and then fuse them into a phase gradient exponential map using the root mean square method: ; In the formula, The phase gradient index; These represent the phase gradient maps after median filtering, representing the four diagonal directions: southeast, northeast, southwest, and northwest. The phase gradient exponential map is normalized to the [0,1] interval to obtain the phase gradient map.

[0008] In some optional embodiments, obtaining the deformation rate map includes the following steps: Obtain the wrapped phase of multiple interferograms, and the phase composition of each interferogram is as follows: ; In the formula, Let x be the winding phase in the i-th interferogram, and let x and y be the row and column number of a pixel in the interferogram, respectively. For terrain phase, For deformation phase, For atmospheric delayed phase, Noise phase; An external digital elevation model (DEM) is introduced for simulation, and the terrain phases are differentially removed. ; Residual terrain errors were removed by constructing a model relating elevation error to spatial baselines, and atmospheric delay phase was simulated and removed by constructing a spatial polynomial model. ; Adaptive filtering is applied to the processed interferogram to suppress the noise phase. Thus, an unwrapped phase diagram dominated by deformation phase is obtained. ; Based on the untangled phase diagram The annual deformation rate of each pixel is calculated using a weighted least squares model, along with the corresponding time intervals. ; In the formula, The deformation rate is in the radar line-of-sight direction; The time interval between the two images used to generate the i-th interferogram is given by n, where n is the number of interferograms. The deformation rate map is generated based on the deformation rate of all pixels.

[0009] In some optional embodiments, the formula for calculating the local adaptive threshold is: ; In the formula, For the phase gradient map at position Pixel value at; The weighting function is Gaussian. The formula for the Gaussian weighting function is: ; In the formula, The standard deviation of the Gaussian function is used to control the strength of the weight distribution and the effective size of the window. Let m be the radius of the neighborhood window, defining the local neighborhood range; m and n are the pixels within the window relative to the center pixel. The offset.

[0010] In some optional embodiments, after classifying the pixels into landslide region pixels and background region pixels, morphological filling processing is performed on the phase discontinuity regions within the landslide region pixels, including the following steps: For the internal regions of landslide distributions that are not effectively identified due to small phase gradient differences caused by uniform deformation, the hole filling algorithm is used to set the pixels classified as background in the landslide region as foreground, thereby generating complete and continuous landslide region pixels. Based on the pixels of the landslide area, multiple geometric bounding boxes are generated, each covering a different landslide area; the corresponding deformation rate distribution map is extracted from the deformation rate map using each geometric bounding box.

[0011] In some optional embodiments, determining the segmentation threshold for each deformation rate distribution map includes the following steps: On the grayscale histogram of each deformation rate distribution map, a straight line is constructed from the position corresponding to the highest peak of the histogram to the brightness value corresponding to the brightest side of the histogram; Starting from the brightness value position corresponding to the brightest side, calculate the vertical distance from each brightness value in the histogram to the straight line; The position of the histogram brightness value corresponding to the maximum vertical distance is the segmentation threshold.

[0012] In some optional embodiments, the deformation region is further subjected to weak deformation enhancement, the specific steps of which are as follows: If the deformation area exhibits a low landslide deformation rate, resulting in grayscale values ​​close to the background, or insufficient local contrast leading to incomplete or blurred boundary segmentation, then an adaptive histogram equalization method with limited contrast is used to process the deformation area.

[0013] In some optional embodiments, the morphological processing and merging of each deformed region to obtain a result map indicating the spatial distribution of the landslide area includes the following steps: Spatial resolution enhancement, morphological combination, and filtering are performed on the deformation area to obtain the vector boundary of the deformation area. After smoothing the boundary and merging, a result map indicating the spatial distribution of the landslide area is obtained. The morphological combination includes dilation, erosion, opening operation, and closing operation.

[0014] Embodiments of the present invention also provide a computer device, comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the landslide identification method described above.

[0015] Embodiments of the present invention also provide a computer-readable storage medium storing a computer program, which, when run by a processor, is capable of executing the landslide identification method described above.

[0016] The landslide identification method provided by this invention has at least the following beneficial effects: The technical solution provided by this invention achieves high-precision automatic landslide identification and boundary localization. It integrates phase gradient maps and deformation rate maps for collaborative analysis. First, the sensitivity of the phase gradient map to deformation boundaries is used for preliminary localization and indication. Then, guided by the geometric bounding boxes generated for each indicated area, the extracted independent deformation rate distribution maps are analyzed in detail. An adaptive threshold segmentation method based on the gray-level histogram features of each distribution map is used to accurately extract significant deformation areas. Finally, morphological processing optimizes the boundaries to obtain a final result that accurately reflects the spatial distribution of landslides. This method effectively combines the boundary indication capability of the phase gradient map with the morphological characterization capability of the deformation rate map, significantly improving the completeness of identification and the accuracy of boundaries. Attached Figure Description

[0017] The accompanying drawings, which are included to provide a further understanding of the invention and form part of this invention, illustrate exemplary embodiments of the invention and are used to explain the invention, but do not constitute an undue limitation of the invention. In the drawings:

[0018] Figure 1 This is a flowchart of a landslide identification method according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the triangulation method provided according to an embodiment of the present invention; Figure 3 This is a deformation rate distribution diagram of Stacking processing according to an embodiment of the present invention; Figure 4 This is an IPGS processing phase gradient indication map provided according to an embodiment of the present invention; Figure 5 This is a phase gradient indication diagram after mean filtering according to an embodiment of the present invention; Figure 6 This is an image showing the result of adaptive threshold binarization processing according to an embodiment of the present invention; Figure 7 This is a distribution map of suspected landslide areas after internal filling treatment according to an embodiment of the present invention; Figure 8 This is a distribution map of the bounding rectangles of each suspected landslide area according to an embodiment of the present invention; Figure 9 This is a schematic diagram of a single landslide deformation rate grayscale mask set according to an embodiment of the present invention; Figure 10 This is a grayscale mask image of the deformation rate of a composite landslide 1 according to an embodiment of the present invention; Figure 11 This is a grayscale mask image of the deformation rate of a non-uniform deformation landslide 2 according to an embodiment of the present invention; Figure 12 This is a gray-scale normalized [0, 255] result diagram of a composite landslide 1 according to an embodiment of the present invention; Figure 13 This is a gray-level normalization [0, 255] result diagram of a non-uniform deformation landslide 2 according to an embodiment of the present invention; Figure 14 This is a graph showing the grayscale normalization [0, 255] result of a non-uniform deformation landslide 2 after weak deformation enhancement according to an embodiment of the present invention. Figure 15 This is a threshold segmentation result diagram of a composite landslide using the triangulation method according to an embodiment of the present invention; Figure 16 This is a threshold segmentation result diagram of a non-uniform deformation landslide using the triangulation method according to an embodiment of the present invention; Figure 17 This is a threshold segmentation result of the triangular method after weak deformation enhancement of a non-uniform deformation landslide 2 according to an embodiment of the present invention; Figure 18 This is a composite landslide threshold segmentation histogram provided according to an embodiment of the present invention; Figure 19 This is a threshold segmentation histogram of a non-uniform deformation landslide using the triangulation method, provided according to an embodiment of the present invention. Figure 20 This is a flowchart of morphological processing and boundary optimization of landslide 1 according to an embodiment of the present invention; Figure 21 This is a flowchart of morphological processing and boundary optimization of landslide 2 according to an embodiment of the present invention; Figure 22 This is a schematic diagram of the phase gradient distribution of landslide hazards according to an embodiment of the present invention; Figure 23 This is a schematic diagram of the deformation rate distribution of landslide hazards according to an embodiment of the present invention; Figure 24 This is a schematic diagram of the flow framework of a landslide identification method according to an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0020] One embodiment of the present invention relates to a landslide identification method. The implementation details of the landslide identification method of this embodiment are described in detail below. The following content is only for the convenience of understanding and is not necessary for implementing this solution.

[0021] The specific process of the landslide identification method in this embodiment can be as follows: Figure 1 As shown, it includes: Step 101: Obtain the phase gradient map and deformation rate map of the target region within a preset time period; Collect Sentinel-1 rising and falling orbit SAR data of the target area, download the corresponding precise orbit data and DEM, and obtain the single-look complex image SLC required for processing through preprocessing. The phase gradient map is obtained through the following steps: A preset number of SAR images are registered to the same master image. Multiple interferometric combinations are generated based on the spatiotemporal baseline threshold. The wrapped phase representation of each interferogram is obtained through interferometric processing: ; In the formula, The winding phase in the i-th interferogram. x and y These represent the number of rows and columns of a pixel in the interferogram, respectively. For terrain phase, For deformation phase, For atmospheric delayed phase, Noise phase; From the interferograms, select M interferograms, and for each interferogram, calculate its phase gradient along at least four directions, where the phase gradient along the southeast direction with step size s is calculated as follows: ; In the formula, For the i-th interferogram at pixel The phase gradient value along the southeast (SE) direction; and For the i-th interferogram in coordinates and The winding phase value at the specified point; s is the step size used when calculating the gradient. The phase gradients of the M interferograms in the same direction are superimposed in the time dimension to obtain the stacked phase gradient in that direction: ; In the formula, For in pixels At this location, the phase gradient value after stacking along the southeast direction; M is the number of high-quality interferograms processed by phase gradient stacking. Median filtering is applied to the stacked phase gradients of all pixels to obtain the filtered phase gradient map: ; In the formula, This is the phase gradient map after median filtering along the southeast direction; (.) represents the median filter function; Repeat the above steps to obtain filtered phase gradient maps in at least four directions, and then fuse them into a phase gradient exponential map using the root mean square method: ; In the formula, The phase gradient index; These represent the phase gradient maps after median filtering, representing the four diagonal directions: southeast, northeast, southwest, and northwest. The phase gradient exponential map is normalized to the [0,1] interval to obtain the phase gradient map.

[0022] Deformation rate maps are obtained through the following steps: Obtain the wrapped phase of multiple interferograms, and the phase composition of each interferogram is as follows: ; In the formula, Let x be the winding phase in the i-th interferogram, and let x and y be the row and column number of a pixel in the interferogram, respectively. For terrain phase, For deformation phase, For atmospheric delayed phase, Noise phase; An external digital elevation model (DEM) is introduced for simulation, and the terrain phases are differentially removed. ; Residual terrain errors were removed by constructing a model relating elevation error to spatial baselines, and atmospheric delay phase was simulated and removed by constructing a spatial polynomial model. ; Adaptive filtering is applied to the processed interferogram to suppress the noise phase. Thus, an unwrapped phase diagram dominated by deformation phase is obtained. ; Based on the untangled phase diagram The annual deformation rate of each pixel is calculated using a weighted least squares model, along with the corresponding time intervals. ; In the formula, The deformation rate is in the radar line-of-sight direction; The time interval between the two images used to generate the i-th interferogram is given by n, where n is the number of interferograms. The deformation rate map is generated based on the deformation rate of all pixels.

[0023] The phase gradient map is subjected to mean filtering to smooth the data, suppress noise, and fill in empty regions.

[0024] Step 102: Based on the features of neighboring pixels and the pixel value of each pixel in the phase gradient map, calculate the corresponding local adaptive threshold for each pixel. By comparing the pixel value of each pixel with the local adaptive threshold of the pixel itself, the pixels are classified into landslide area pixels and background area pixels. The formula for calculating the local adaptive threshold is: ; In the formula, For the phase gradient map at position Pixel value at; The weighting function is Gaussian. The formula for the Gaussian weighting function is: ; In the formula, The standard deviation of the Gaussian function is used to control the strength of the weight distribution and the effective size of the window. Let m be the radius of the neighborhood window, defining the local neighborhood range; m and n are the pixels within the window relative to the center pixel. The offset.

[0025] Morphological internal filling processing is performed on phase discontinuity regions within the landslide area pixels. This morphological internal filling processing targets the internal regions of the landslide distribution where the phase gradient difference is small due to uniform deformation and they are not effectively identified. By using a hole filling algorithm, pixels classified as background in these regions are set as foreground pixels, thereby generating complete and continuous landslide area pixels. Based on these landslide area pixels, multiple geometric bounding boxes are generated, each covering a different landslide area. The corresponding deformation rate distribution map is extracted from the deformation rate map using each geometric bounding box.

[0026] Step 103: Based on the pixels of the landslide area, generate multiple geometric bounding boxes that cover different landslide areas respectively; Step 104: Extract the corresponding deformation rate distribution map from the deformation rate map using each geometric bounding box; The deformation rate distribution map is taken as an absolute value, and the grayscale data after taking the absolute value is mapped to the interval [0, 255] using the MinMax normalization method; where, for a grayscale value X of a certain point in the grayscale map after taking the absolute value, the normalized value is... Calculated using the following formula: ; In the formula, and These are the maximum and minimum grayscale values ​​in the grayscale image after taking their absolute values.

[0027] Step 105: Based on the grayscale histogram features of each deformation rate distribution map, determine the segmentation threshold of each deformation rate distribution map; segment the corresponding deformation rate distribution map according to the segmentation threshold, and extract the deformation regions in the deformation rate distribution map that are higher than the segmentation threshold. The segmentation threshold is determined as follows: On the grayscale histogram of each deformation rate distribution map, a straight line is constructed from the position corresponding to the highest peak of the histogram to the brightness value corresponding to the brightest side of the histogram; starting from the brightness value position corresponding to the brightest side, the vertical distance from each brightness value in the histogram to the straight line is calculated; the brightness value position of the histogram corresponding to the maximum value of the vertical distance is the segmentation threshold.

[0028] A diagram of the trigonometric method is shown below. Figure 2 As shown in the figure, a straight line is constructed on the histogram from the highest peak to the brightest corresponding histogram brightness value Bmax. Starting from Bmax, the vertical distance from each corresponding histogram B to the straight line is calculated until d is maximized. The histogram position corresponding to the maximum distance is the threshold corresponding to image binarization.

[0029] If the deformation area exhibits a low landslide deformation rate, resulting in grayscale values ​​close to the background, or insufficient local contrast leading to incomplete or blurred boundary segmentation, then an adaptive histogram equalization method with limited contrast is used to process the deformation area.

[0030] The steps of the contrast-limited adaptive histogram equalization method are as follows: The pixels in the deformed region are evenly divided into their respective gray levels, and the average number of pixels allocated to each gray level is calculated. The calculation formula is: ; In the formula, This represents the number of pixels in the image along the x-axis. This represents the number of pixels in the image along the y-axis. This represents the number of gray levels contained in the image. The threshold is then determined. for:

[0031] ; In the formula, To truncate to the maximum possible value, the final pixel count exceeded the threshold. The portion will be redistributed:

[0032] ; in, The total number of pixels captured. This assigns an average number of pixels to each grayscale level. This adaptive adjustment process stretches the dynamic range of the image, enhances image contrast, and highlights weak areas of surface deformation. Finally, histogram equalization is performed on each image after pixel redistribution to achieve the effect of enhancing weakly deformed areas.

[0033] Step 106: Morphological processing and merging of each deformed region are performed to obtain a result map indicating the spatial distribution of the landslide area.

[0034] Spatial resolution enhancement, morphological combination, and filtering are performed on the deformation area to obtain the vector boundary of the deformation area. After smoothing the boundary and merging, a result map indicating the spatial distribution of the landslide area is obtained. The morphological combination includes dilation, erosion, opening operation, and closing operation.

[0035] Morphological image processing involves moving a structuring element (rectangle, circle, ellipse, cross, etc.) within an image and then performing set operations such as intersection and union with the underlying binary image. By defining and selecting the structuring element for comparison and matching, the shape and structure of the image to be processed are obtained. Then, operations from digital morphology, such as dilation, erosion, opening, and closing, are introduced. Through a combination of morphological operators, background interference is removed, the target geometry is determined, and jagged edges and internal null values ​​in the slope are resolved, thereby ensuring the spatial structural features of the image and the target geometry.

[0036] The dilation operation expands the image boundaries, filling in internal holes, and typically makes the image brighter after processing. The grayscale dilation operation is defined as follows:

[0037] ; In the formula It is a grayscale image. For structuring elements that require defining an origin, , They are , The domain of the pixel In the image When it moves up, it will move the image The maximum value on is assigned to The origin is used to eliminate noise in the feature region, fill in null values, and smooth lines.

[0038] Erosion operations shrink image boundaries, removing isolated pixels and typically resulting in a darker image. The grayscale erosion operation is defined as follows:

[0039] ; when exist When moving upwards, Included The minimum value of the pixel is assigned to This removes background noise and retains only the effective features.

[0040] The opening operation first performs erosion on the image, followed by dilation. This step removes isolated pixels, connections between pixels, and rough edges, smoothing the image contours. Its calculation is defined as follows:

[0041] ; Gray-level closing and gray-level opening are dual operations. Closing first performs dilation on the image, then erosion. This step smooths the image, but unlike opening, it fills in null values ​​and increases the spatial correlation of pixels. Its calculation is defined as follows:

[0042] ; The processing also includes enhancing the image spatial resolution and introducing median filtering to remove existing noise and introduced interference. After obtaining the landslide identification result through combined operations, if a simplified surface is not selected during the raster-to-surface conversion process, the boundary will be completely consistent with the pixel edge. Therefore, after removing individual deformation holes, the vector surface needs to be smoothed. Finally, the identification result is verified against the optical image to confirm it as a potential landslide hazard.

[0043] Example 1: This patent takes a landslide-prone area in a typical watershed as the research area, and selects a composite landslide (marked as landslide 1) and a non-uniform deformation landslide (marked as landslide 2) in the area as specific embodiments, and applies the described technical solution to carry out an automatic identification process.

[0044] For large-scale landslide identification, Improved Phase Gradient Stacking (IPGS) is used to suppress noise interference and avoid errors that may be introduced by phase unwrapping. Meanwhile, the InSAR Stacking method can more accurately reflect the true deformation range and morphology of the landslide by quantifying the deformation rate. To ensure the accuracy of landslide identification and minimize error interference, this implementation adopts a combined strategy: first, the distribution location indication of the landslide hazard is obtained based on the IPGS processing results; then, under the indication of this spatial location, the high-precision deformation rate results obtained by Stacking technology are used for fine segmentation, thereby ultimately determining the accurate boundary range of the landslide.

[0045] Specifically, the execution of this embodiment corresponds to steps one through four in the technical solution. Its core objective is to extract a grayscale deformation rate map (i.e., a deformation rate grayscale mask) of the corresponding region, based on the deformation anomaly area identified by IPGS as a spatial guide, thereby completing the boundary extraction of a single landslide hazard. The processing results and process correspondence of typical landslide hazard areas are shown in Table 1. This table clearly shows the key intermediate results from the initial IPGS indication to the acquisition of the deformation rate mask.

[0046] Table 1. Correspondence between the treatment results and procedures for landslide hazard areas Two typical deformation rate distribution maps (composite landslide 1 and non-uniform deformation landslide 2) are selected to illustrate subsequent processing. Among them, weak deformation enhancement was performed on non-uniform landslide 2, resulting in a more complete segmentation compared to the unenhanced version. The deformation rate distribution map segmentation is shown in Table 2:

[0047] Table 2 Deformation rate distribution map segmentation Morphological boundary optimization stage: When optimizing the boundaries of the segmentation results, the image spatial resolution is first enhanced, and digital morphological operations such as dilation, erosion, opening, and closing are introduced to smooth and optimize the boundaries. Isolated pixels are removed through filtering. This example provides a morphological processing scheme for complex and non-uniform landslides, which can be adjusted according to the actual effect in large-scale implementation.

[0048] The flowchart of morphological processing and boundary optimization for landslide 1 is as follows: Figure 20 As shown, Figure 20 Taking landslide 1 as an example, the process of optimizing the boundary of a landslide area using morphological operations is demonstrated: First, an opening operation is performed on the suspected landslide area obtained from the initial segmentation to remove minor noise; then, isolated pixels are removed through filtering; finally, the boundary is optimized through morphological smoothing, resulting in a landslide area that closely matches the actual shape (IoU reaches 0.9190). This process is suitable for the refinement of the boundaries of both complex and uniform landslides, and the morphological operation parameters can be adjusted according to the actual scenario.

[0049] The flowchart of morphological processing and boundary optimization for landslide 2 is as follows: Figure 21 As shown, Figure 21Taking landslide 2 as an example, this paper demonstrates a morphological processing workflow suitable for non-uniform landslide areas with dispersed noise: First, a closing operation is performed on the initially segmented suspected landslide areas to fill the small cavities inside the areas; then, filtering is used to remove peripheral dispersed noise points; subsequently, a dilation operation is performed to integrate the associated areas; finally, boundary smoothing is performed to obtain a precise landslide area boundary with an IoU of 0.8994. This workflow is suitable for non-uniform landslide types with internal cavities and peripheral dispersed noise, and the operating parameters can be adjusted according to the actual landslide characteristics.

[0050] Obtain the phase gradient map and deformation rate map of the target area within a preset time period. A schematic diagram of the phase gradient distribution of landslide hazard is shown below. Figure 22 As shown in the figure, this is a visualization of the distribution of landslide hazards in a certain area. Figure 22 The image shows the visualization results of landslide hazard distribution in a certain area. The black solid line in the image represents the boundary of the landslide hazard area after being identified and optimized by this method. Combined with geographic coordinates and scale (0-5 km), the spatial location and distribution range of each landslide hazard can be presented intuitively.

[0051] A schematic diagram of the deformation rate distribution of landslide hazards is shown below. Figure 23 As shown in the figure, the identified landslide hazard area (black solid line) is overlaid with the original remote sensing image of the area. This allows for a clear comparison of the location and surface features of the landslide hazard area in the actual geographical scene, verifying the matching of the identification results with the actual terrain, and facilitating subsequent on-site verification and risk assessment.

[0052] To verify the effectiveness and reliability of the automatic landslide identification method described in this patent, a composite landslide 1 and a non-uniform deformation landslide 2 within the study area were selected as typical examples for experimental verification. The implementation process strictly followed the steps of the described technical solution. First, the Sentinel-1 time-series SAR data was processed using the IPGS method to generate a phase gradient map. This method effectively suppresses systematic errors such as atmospheric delay and orbital residuals, as well as random noise, by calculating and superimposing the phase gradients of multi-temporal interferograms. This allows for robust identification of the initial spatial location of potential landslide hazards, avoiding misjudgments that may be caused by unwrapping errors in traditional InSAR.

[0053] Based on the spatial indication of the IPGS results, the annual average deformation rate map, which reflects the actual surface movement trend and is obtained using stacking technology, was used to independently analyze each suspected landslide area. Through triangulation threshold segmentation, areas of continuous deformation were extracted from the deformation rate map, thereby obtaining the true shape and distribution range of landslide hazards in the deformation dimension.

[0054] To further optimize the identification boundary, digital morphological processing, including dilation, erosion, and opening / closing operations, was introduced into the segmented binary deformed region. This step effectively removed isolated noise pixels from the result, smoothed irregular edges, and filled internal voids caused by uneven deformation, ultimately yielding an optimized result that highly matches the boundary of the natural landslide landform. Quantitative evaluation showed that the intersection-over-union (IoU) of the landslide region boundary after the above processing was significantly improved compared to the ground truth. For the example landslides, the IoU values ​​all reached over 85% (e.g., the IoU of landslide 1 was 0.9190), demonstrating the effectiveness of the boundary optimization.

[0055] In summary, this embodiment verifies the automatic landslide identification method integrating IPGS, Stacking, and morphological processing. This method eliminates the need for manual sample annotation or complex machine learning model training, directly utilizing interferometric phase information to achieve a fully automated process from hazard location and deformation zone extraction to boundary optimization. It features high processing efficiency and wide applicability, making it suitable for systematic landslide hazard surveys and early identification in large-scale areas, providing a reliable technical tool for regional landslide hazard assessment and early warning.

[0056] The framework of landslide identification method process is as follows: Figure 24 As shown in the figure, this diagram presents the complete architecture of the wide-area landslide InSAR automatic identification method based on digital morphology. It consists of four core modules: First, SAR data processing. This module takes ascending SAR imagery, SRTM DEM, descending SAR imagery, and precise orbit data as input, generating phase gradient maps through IPGS processing and deformation rate maps through Stacking InSAR processing. Second, deformation region extraction. This module sequentially performs mean filtering, adaptive threshold binarization, and internal filling operations on the phase gradient map to generate an outer bounding box, which is then combined with the deformation rate map to complete a grayscale mask. Third, weak deformation enhancement. This module first performs grayscale normalization on the data, then extracts the deformation region through triangulation threshold segmentation. If local contrast is insufficient, the CLAHE method is used for enhancement. Fourth, morphological processing. This module first enhances the resolution of the data, then iteratively performs morphological operations such as dilation, erosion, opening, closing, and filtering. The results are then subjected to hole filling, boundary smoothing, and merging. Finally, the results are combined with optical imagery for verification, outputting accurate landslide identification results.

[0057] The steps of the various methods described above are only for clarity. In practice, they can be combined into one step or some steps can be split into multiple steps. As long as they include the same logical relationship, they are all within the protection scope of this invention. Adding insignificant modifications or introducing insignificant designs to the algorithm or process, without changing the core design of the algorithm and process, are also within the protection scope of this invention.

[0058] Another embodiment of the present invention relates to a computer-readable storage medium storing a computer program. When executed by a processor, the computer program implements the method embodiments described above.

[0059] That is, those skilled in the art will understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware. This program is stored in a storage medium and includes several instructions to cause a device (which may be a microcontroller, chip, etc.) or processor to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0060] Those skilled in the art will understand that the above embodiments are specific embodiments for implementing the present invention, and in practical applications, various changes can be made to them in form and detail without departing from the spirit and scope of the present invention.

Claims

1. A landslide identification method, characterized in that, The method includes: Acquire the phase gradient map and deformation rate map of the target region within a preset time period; Based on the features of neighboring pixels and the pixel value of each pixel in the phase gradient map, a corresponding local adaptive threshold is calculated for each pixel. By comparing the pixel value of each pixel with the local adaptive threshold of the pixel itself, the pixels are classified into landslide area pixels and background area pixels. Based on the pixels of the landslide area, generate multiple geometric bounding boxes that cover different landslide areas; The corresponding deformation rate distribution map is extracted from the deformation rate map using each geometric bounding box; Based on the grayscale histogram features of each deformation rate distribution map, the segmentation threshold of each deformation rate distribution map is determined; the corresponding deformation rate distribution map is segmented according to the segmentation threshold, and the deformation regions in the deformation rate distribution map that are higher than the segmentation threshold are extracted. Morphological processing and merging of each deformation region yielded a map showing the spatial distribution of the landslide area.

2. The landslide identification method as described in claim 1, characterized in that, The acquisition of the phase gradient map includes the following steps: A preset number of SAR images are registered to the same master image. Multiple interferometric combinations are generated based on the spatiotemporal baseline threshold. The wrapped phase representation of each interferogram is obtained through interferometric processing: ； In the formula, The winding phase in the i-th interferogram. x and y These represent the number of rows and columns of a pixel in the interferogram, respectively. For terrain phase, For deformation phase, For atmospheric delayed phase, Noise phase; From the interferograms, select M interferograms, and for each interferogram, calculate its phase gradient along at least four directions, where the phase gradient along the southeast direction with step size s is calculated as follows: ； In the formula, For the i-th interferogram at pixel The phase gradient value along the southeast (SE) direction; and For the i-th interferogram in coordinates and The winding phase value at the specified point; s is the step size used when calculating the gradient. The phase gradients of the M interferograms in the same direction are superimposed in the time dimension to obtain the stacked phase gradient in that direction: ； In the formula, For in pixels At this location, the phase gradient value after stacking along the southeast direction; M is the number of high-quality interferograms processed by phase gradient stacking. Median filtering is applied to the stacked phase gradients of all pixels to obtain the filtered phase gradient map: ； In the formula, This is the phase gradient map after median filtering along the southeast direction; (.) represents the median filter function; Repeat the above steps to obtain filtered phase gradient maps in at least four directions, and then fuse them into a phase gradient exponential map using the root mean square method: ； In the formula, The phase gradient index; These represent the phase gradient maps after median filtering, representing the four diagonal directions: southeast, northeast, southwest, and northwest. The phase gradient exponential map is normalized to the [0,1] interval to obtain the phase gradient map.

3. The landslide identification method as described in claim 1, characterized in that, The acquisition of the deformation rate map includes the following steps: Obtain the wrapped phase of multiple interferograms, and the phase composition of each interferogram is as follows: ； In the formula, Let x be the winding phase in the i-th interferogram, and let x and y be the row and column number of a pixel in the interferogram, respectively. For terrain phase, For deformation phase, For atmospheric delayed phase, Noise phase; An external digital elevation model (DEM) is introduced for simulation, and the terrain phases are differentially removed. ; Residual terrain errors were removed by constructing a model relating elevation error to spatial baselines, and atmospheric delay phase was simulated and removed by constructing a spatial polynomial model. ; Adaptive filtering is applied to the processed interferogram to suppress the noise phase. Thus, an unwrapped phase diagram dominated by deformation phase is obtained. ; Based on the untangled phase diagram The annual deformation rate of each pixel is calculated using a weighted least squares model, along with the corresponding time intervals. ； In the formula, The deformation rate is in the radar line-of-sight direction; The time interval between the two images used to generate the i-th interferogram is given by n, where n is the number of interferograms. The deformation rate map is generated based on the deformation rate of all pixels.

4. The landslide identification method as described in claim 1, characterized in that, The formula for calculating the local adaptive threshold is: ； In the formula, For the phase gradient map at position Pixel value at; The weighting function is Gaussian. The formula for the Gaussian weighting function is: ； In the formula, The standard deviation of the Gaussian function is used to control the strength of the weight distribution and the effective size of the window. Let m be the radius of the neighborhood window, defining the local neighborhood range; m and n are the pixels within the window relative to the center pixel. The offset.

5. The landslide identification method as described in claim 1, characterized in that, The method also includes, after classifying pixels into landslide region pixels and background region pixels, performing morphological internal filling processing on phase discontinuity regions within landslide region pixels, including the following steps: For the internal regions of landslide distributions that are not effectively identified due to small phase gradient differences caused by uniform deformation, the hole filling algorithm is used to set the pixels classified as background in the landslide region as foreground, thereby generating complete and continuous landslide region pixels. Based on the pixels of the landslide area, multiple geometric bounding boxes are generated, each covering a different landslide area; the corresponding deformation rate distribution map is extracted from the deformation rate map using each geometric bounding box.

6. The landslide identification method as described in claim 1, characterized in that, Determining the segmentation threshold for each deformation rate distribution map includes the following steps: On the grayscale histogram of each deformation rate distribution map, a straight line is constructed from the position corresponding to the highest peak of the histogram to the brightness value corresponding to the brightest side of the histogram; Starting from the brightness value position corresponding to the brightest side, calculate the vertical distance from each brightness value in the histogram to the straight line; The position of the histogram brightness value corresponding to the maximum vertical distance is the segmentation threshold.

7. The landslide identification method as described in claim 1, characterized in that, It also includes weak deformation enhancement of the deformed region, including the following steps: If the deformation area exhibits a low landslide deformation rate, resulting in grayscale values ​​close to the background, or insufficient local contrast leading to incomplete or blurred boundary segmentation, then an adaptive histogram equalization method with limited contrast is used to process the deformation area.

8. The landslide identification method as described in claim 1, characterized in that, The morphological processing and merging of each deformation region to obtain a map indicating the spatial distribution of the landslide area includes the following steps: Spatial resolution enhancement, morphological combination, and filtering are performed on the deformation area to obtain the vector boundary of the deformation area. After smoothing the boundary and merging, a result map indicating the spatial distribution of the landslide area is obtained. The morphological combination includes dilation, erosion, opening operation, and closing operation.

9. A computer system, characterized in that, include: At least one processor; And a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the landslide identification method as described in any one of claims 1 to 8.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, is capable of performing the landslide identification method as defined in any one of claims 1 to 8.

Images