An InSAR technology-based landslide partition automatic processing method

Abstract

The application provides a landslide partition automatic processing method based on InSAR technology, belongs to the technical field of landslide deformation remote sensing monitoring, and aims at the problems that the existing single InSAR method is difficult to adapt to multiple ground objects in a complex landslide area, and is prone to uneven monitoring point density, missing deformation in a low-coherence area and inconsistent partition results. Synthetic aperture radar images, orbit ephemeris, digital elevation models, land cover vector data and landslide geological survey data are acquired, and the coherence quality index, the intensity quality index, the radar visibility index and the spatiotemporal baseline suitability are calculated. The land cover type and the landslide geological partition are superimposed to demarcate a monitoring subarea, the processing adaptation index is formed in combination with the scatterer distribution characteristics, and the InSAR processing strategy is matched. The partition timing deformation inversion is carried out based on a unified reference datum, and the timing deformation monitoring result of the whole landslide area is obtained after processing. The application improves the integrity, consistency and automation level of deformation monitoring in a complex landslide area.

Description

Technical Field

[0001] This invention belongs to the field of landslide deformation remote sensing monitoring technology, specifically relating to a method for automatic landslide zoning using synthetic aperture radar interferometry. Background Technology

[0002] Landslides are a common type of geological hazard in mountainous, hilly, and engineering-disturbed areas. Their occurrence and development are typically influenced by a combination of factors, including topography, soil and rock structure, rainfall infiltration, river erosion, engineering excavation, and changes in vegetation cover. Before instability, landslides often undergo stages of slow deformation, localized accelerated deformation, and overall instability. Therefore, continuous, detailed, and large-scale surface deformation monitoring of landslide areas is a crucial technical aspect of landslide hazard identification, disaster early warning, and risk prevention. Traditional landslide monitoring methods typically include total stations, Global Navigation Satellite Systems (GNSS), fissure gauges, inclinometers, and manual inspections. While these methods can obtain high-precision observation data at local monitoring points, their monitoring range is limited, deployment costs are high, and they are easily constrained by terrain accessibility, monitoring point density, and on-site maintenance conditions. Consequently, they cannot meet the needs of large-scale, long-term, and high-frequency monitoring of landslides in complex mountainous areas.

[0003] Interferometric Synthetic Aperture Radar (InSAR) technology can acquire surface deformation information using multi-phase Synthetic Aperture Radar (SAR) imagery. It boasts advantages such as wide coverage, strong all-weather and all-time observation capabilities, high spatial resolution, and suitability for long-term sequence analysis, and has been gradually applied in the field of landslide deformation monitoring. With the development of time-series InSAR technology, methods such as Persistent Scatterer Interferometric Synthetic Aperture Radar (PS-InSAR), Distributed Scatterer Interferometric Synthetic Aperture Radar (DS-InSAR), and Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) have been used to extract deformation in regions with different scattering characteristics. Among them, PS-InSAR is suitable for areas with many stable strong scattering targets such as artificial buildings and exposed rock masses; DS-InSAR is suitable for areas with a certain number of homogeneous distributed scatterers; SBAS-InSAR is suitable for reducing the decoherence effect by combining short-time and spatial baselines and is used for time-series deformation inversion under complex land cover conditions.

[0004] However, in actual landslide monitoring scenarios, landslide areas often simultaneously contain multiple land cover types, such as exposed bedrock, loose deposits, vegetated areas, roads, and buildings. The radar backscattering characteristics, coherence stability, and distribution of identifiable scatterers vary significantly across different areas. Existing processing methods often employ a single InSAR method to uniformly process the entire landslide area, or rely on manual experience to select processing methods for different areas. This easily leads to mismatched processing strategies for high-coherence and low-coherence areas, resulting in uneven monitoring point density, insufficient effective monitoring points in landslide deposit areas, and missing deformation information in vegetated areas. Furthermore, if deformation inversion is performed independently in different zones, the reference points, deformation benchmarks, and error levels between zones may be inconsistent, easily causing abrupt deformation changes or systematic biases at zone boundaries, affecting the continuity and reliability of deformation results for the entire landslide area. Furthermore, existing technologies are insufficient in their joint utilization of land cover data, landslide geological zoning, coherence indices, intensity indices, and scatterer distribution characteristics. They lack a complete processing workflow that extends from automatic delineation of monitoring sub-regions, automatic matching of InSAR processing strategies, construction of a unified reference benchmark, to error correction and spatial fusion of zoning results.

[0005] To address this, an automatic InSAR landslide zoning processing method is proposed that can adapt to complex land cover and landslide geological zoning characteristics, thereby improving the completeness, consistency, and automation of landslide temporal deformation monitoring results. Summary of the Invention

[0006] The present invention aims to solve at least one of the technical problems existing in the prior art, and to provide an automatic landslide zoning processing method based on InSAR technology.

[0007] This invention provides an automatic landslide zoning processing method based on InSAR technology, comprising the following steps: S1: Acquire SAR images, orbital ephemeris data, digital elevation model data, land cover vector data, and landslide geological survey data for the monitoring period in the study area; S2: Based on the orbital ephemeris data and the digital elevation model data, the SAR image is preprocessed, and the coherence quality index and intensity quality index for zoning processing are calculated; the land cover vector data is spatially overlaid with the landslide geological zoning in the landslide geological survey data to automatically delineate monitoring sub-zones; a radar visibility mask is generated based on the observation geometry of the digital elevation model data and the SAR image, and the radar visibility index is determined according to the distribution of the radar visibility mask in each monitoring sub-zone; and the spatiotemporal baseline suitability is determined according to the temporal and spatial baseline distribution of the SAR image. S3: Based on the land cover type, coherence quality index, intensity quality index, scatterer distribution characteristics, radar visibility index and spatiotemporal baseline suitability of each monitoring sub-region, a processing adaptation index is formed for each monitoring sub-region, and the corresponding InSAR processing strategy and scatterer screening rules are matched based on the processing adaptation index. S4: Based on the InSAR processing strategy and scatterer screening rules obtained by matching each monitoring sub-region, and combined with the deformation reference benchmark established by the unified reference point or reference point set, perform time-series deformation inversion on each monitoring sub-region to obtain the time-series deformation results of each monitoring sub-region. S5: Spatial registration, system error correction, outlier removal and spatial seamless fusion are performed on the temporal deformation results of each monitoring sub-region to obtain the temporal deformation monitoring results under a unified reference benchmark for the entire landslide area. The temporal deformation monitoring results include the radar line-of-sight temporal deformation.

[0008] Further, in step S2, the preprocessing includes track correction, image registration, differential interferometry, filtering, and quality index calculation, wherein the quality index calculation includes calculating the average coherence coefficient and average intensity value of the pixels.

[0009] Specifically, in step S2, the spatial overlay includes: selecting a reference SAR image from the SAR image, converting the land cover vector data to the coordinate system corresponding to the reference SAR image, reclassifying the land cover vector data, and overlaying the reclassification result with the landslide geological zone to obtain the boundary of the monitoring sub-zone.

[0010] Specifically, the monitoring sub-areas include exposed bedrock areas, landslide accumulation areas, vegetation-covered areas, and artificial building areas, and are delineated according to the superposition relationship between landslide geological attributes and land cover types.

[0011] Preferably, in step S3, the scatterer distribution characteristics include permanent scatterer candidate point distribution characteristics and distributed scatterer candidate point distribution characteristics. The permanent scatterer candidate point distribution characteristics are determined based on the pixel distribution that satisfies preset coherence and intensity conditions, and the distributed scatterer candidate point distribution characteristics are determined based on the statistical consistency between candidate pixels and neighboring pixels. The matching of the InSAR processing strategy includes: when the permanent scatterer candidate point distribution in the exposed bedrock area or the artificial building area satisfies the deformation solution network construction conditions, the PS-InSAR processing strategy is matched; when the permanent scatterer candidate points and distributed scatterer candidate points in the landslide deposit area satisfy the joint solution conditions, the PS-InSAR combined with DS-InSAR processing strategy is matched; when the vegetation cover area satisfies the short-spatial-temporal baseline differential interferometry pair construction conditions, the SBAS-InSAR processing strategy is matched.

[0012] Specifically, when performing time-series deformation inversion on the exposed bedrock area and the artificial building area, permanent scatterer points are selected based on the average coherence coefficient in the coherence quality index and the average intensity value in the intensity quality index, and a PS-InSAR deformation solution network is constructed based on the permanent scatterer points.

[0013] Furthermore, when performing time-series deformation inversion on the landslide accumulation area, permanent scatterer points are first screened within the landslide accumulation area, and a deformation calculation benchmark is established based on the permanent scatterer points within the landslide accumulation area. Then, distributed scatterer points are screened according to the distribution characteristics of homogeneous scatterers, and joint deformation calculation is performed on the permanent scatterer points and the distributed scatterer points.

[0014] Furthermore, when performing temporal deformation inversion on the vegetation cover area, a short-spatial-temporal baseline differential interferometric pair combination is constructed based on the spatiotemporal baseline constraints, and SBAS-InSAR temporal deformation inversion is performed based on the short-spatial-temporal baseline differential interferometric pair combination.

[0015] Furthermore, in step S4, the unified reference point or reference point set is selected from the overlapping area of ​​the buffer zones of adjacent monitoring sub-regions, and the unified reference point or reference point set satisfies the preset coherence condition, the preset intensity condition, and the preset temporal deformation stability condition.

[0016] Specifically, in step S5, the seamless spatial fusion includes: calculating the partition deformation residual based on the same-named points, adjacent points or grid cells in the overlapping areas of adjacent monitoring sub-regions, performing system error correction and regional network adjustment based on the partition deformation residual, and sampling the corrected temporal deformation results to a unified spatial resolution grid.

[0017] The beneficial effects of this invention are as follows: This invention automatically delineates exposed bedrock areas, landslide deposit areas, vegetation-covered areas, and artificial structure areas by spatially overlaying land cover vector data and landslide geological survey data, enabling the monitoring sub-regions to simultaneously reflect the scattering characteristics of ground features and the geological attributes of landslides. Based on the ground feature type, coherence quality index, intensity quality index, scatterer distribution characteristics, radar visibility index, and spatiotemporal baseline suitability of each monitoring sub-region, this invention automatically matches the corresponding InSAR processing strategy, avoiding manual selection of processing techniques and improving the automation level of the processing flow. For different monitoring sub-regions, this invention employs PS-InSAR, PS-InSAR combined with DS-InSAR, and SBAS-InSAR for temporal deformation inversion, increasing the effective monitoring point density in complex landslide areas and reducing the loss of deformation information in low-coherence areas and loose deposit areas. By establishing a unified deformation reference benchmark through a unified reference point or reference point set, and performing systematic error correction, outlier removal, regional network adjustment, and seamless spatial fusion on the zoning results, the consistency and continuity of deformation results across different zones are improved. Attached Figure Description

[0018] Figure 1 This is a flowchart illustrating the steps of an automatic landslide zoning processing method based on InSAR technology, according to a specific embodiment of the present invention. Figure 2 Spatial distribution map of monitoring points processed using traditional unified PS-InSAR; Figure 3 This is a spatial distribution map of the monitoring points for the automatic partitioning processing of this invention; Figure 4 A comparison chart of the number of effective monitoring points in each monitoring sub-region between traditional unified PS-InSAR processing and the automatic partitioning processing of this invention. Figure 5 A comparison chart of the effective area coverage of each monitoring sub-region between traditional unified PS-InSAR processing and the automatic zoning processing of this invention. Figure 6 A comparison chart of the number of effective monitoring points and the proportion of void area in low coherence areas between traditional unified PS-InSAR processing and the automatic partitioning processing of this invention; Figure 7 A comparison chart of inter-regional deformation residuals and root mean square errors of GNSS verification points between traditional unified PS-InSAR processing and the automatic partitioning processing of this invention; Figure 8 This is a comparison chart of the residual deformation of the partition boundary, the root mean square error of deformation of the overlapping area, and the standard deviation of deformation of the overlapping area before and after partition fusion according to the method of the present invention. Figure 9 A comparison of the number of pixels with boundary jumps before and after partition fusion. Detailed Implementation

[0019] To enable those skilled in the art to better understand the technical solution of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0020] like Figure 1 As shown in the figure, a landslide zoning automatic processing method based on InSAR technology provided by a specific embodiment of the present invention includes the following steps: S1: Acquire SAR images, orbital ephemeris data, digital elevation model data, land cover vector data, and landslide geological survey data for the monitoring period in the study area; S2: Based on the orbital ephemeris data and the digital elevation model data, the SAR image is preprocessed, and the coherence quality index and intensity quality index for zoning processing are calculated; the land cover vector data is spatially overlaid with the landslide geological zoning in the landslide geological survey data to automatically delineate monitoring sub-zones; a radar visibility mask is generated based on the observation geometry of the digital elevation model data and the SAR image, and the radar visibility index is determined according to the distribution of the radar visibility mask in each monitoring sub-zone; and the spatiotemporal baseline suitability is determined according to the temporal and spatial baseline distribution of the SAR image. S3: Based on the land cover type, coherence quality index, intensity quality index, scatterer distribution characteristics, radar visibility index and spatiotemporal baseline suitability of each monitoring sub-region, a processing adaptation index is formed for each monitoring sub-region, and the corresponding InSAR processing strategy and scatterer screening rules are matched based on the processing adaptation index. S4: Based on the InSAR processing strategy and scatterer screening rules obtained by matching each monitoring sub-region, and combined with the deformation reference benchmark established by the unified reference point or reference point set, perform time-series deformation inversion on each monitoring sub-region to obtain the time-series deformation results of each monitoring sub-region. S5: Spatial registration, systematic error correction, outlier removal, and seamless spatial fusion are performed on the temporal deformation results of each monitoring sub-region to obtain the temporal deformation monitoring results under a unified reference benchmark for the entire landslide area.

[0021] In one implementation, PS-InSAR refers to Permanent Scatterer Synthetic Aperture Radar Interferometry; DS-InSAR refers to Distributed Scatterer Synthetic Aperture Radar Interferometry; SBAS-InSAR refers to Small Baseline Assembled Synthetic Aperture Radar Interferometry; the SAR imagery consists of multiple synthetic aperture radar images covering the same study area and within the same monitoring period; the orbital ephemeris data is used to correct orbital errors in the SAR images; the digital elevation model data is used to simulate and remove topographic phase; the land cover vector data is used to identify the land cover type of the study area; and the landslide geological survey data is used to determine the landslide boundary, landslide backwall, landslide deposit, landslide impact zone, and relatively stable area.

[0022] Preferably, the coherence quality index is used to characterize the phase stability of a pixel or candidate scatterer during the monitoring period, and the intensity quality index is used to characterize the radar backscattering intensity and stability of a pixel or candidate scatterer. The coherence quality index and the intensity quality index are used together for subsequent stable point screening, scatterer screening, processing strategy matching, and selection of a unified reference point or reference point set.

[0023] Furthermore, the orbit correction, registration, and differential interferometry processing are performed based on a reference SAR image; the reference SAR image is determined according to the temporal and spatial baseline distribution among multiple SAR images, so that it can serve as a common reference image for subsequent image registration and interferometry processing.

[0024] Based on the above basic implementation method, in step S2, the preprocessing includes orbit correction, image registration, differential interferometry, filtering, and quality index calculation. The orbital ephemeris data includes precise orbit determination ephemeris data for orbit correction or orbital parameter data that can provide satellite orbital position information; the digital elevation model data is used to simulate terrain phase and removes terrain phase during differential interferometry processing; the image registration includes registering other SAR images (excluding the reference SAR image) within the monitoring period to the reference SAR image; the differential interferometry includes forming interferometric pairs based on the registered SAR images and removing terrain phase using the digital elevation model data to obtain a differential interferogram; the filtering is used to suppress phase noise in the differential interferogram. The quality index calculation includes calculating the average coherence coefficient and average intensity value of the pixels. The average coherence coefficient is obtained statistically based on the coherence coefficients of the same pixel in multiple interferometric pairs within the monitoring period, and the average intensity value is obtained statistically based on the backscattering intensity of the same pixel in multiple SAR images. Based on the digital elevation model data, the incident angle and orbital direction of the SAR image, and the slope aspect and gradient of the study area, a radar visibility mask is generated to identify overlapping areas, shadow areas, perspective contraction areas, and low observation sensitivity areas, which are then designated as low-confidence areas. Radar visibility indices are determined based on the area proportion and spatial distribution of these low-confidence areas within each monitoring sub-region of the radar visibility mask. The suitability of the spatiotemporal baseline is determined based on the temporal baseline, spatial baseline, and the number of differential interferometric pairs satisfying the preset spatiotemporal baseline constraints of the SAR image.

[0025] In this embodiment, the radar visibility index and spatiotemporal baseline suitability The method for determining the value is as follows: Radar visibility indicators: ,in The proportion of low-confidence regions. Spatial clustering degree (0~1). For "high", It is rated as "medium". The rating is "low". Spatiotemporal baseline suitability: , To meet d and Number of interference pairs of m. For "high", It is rated as "medium". It is "low"; among them, The sequence number of the monitored sub-region, This serves as the time baseline between two SAR images. The vertical spatial baseline between two SAR images; The number of interferometric pairs required to simultaneously satisfy both temporal and spatial baseline constraints and generate effective differential interferometric results within the i-th monitoring sub-region; 30 represents the preset minimum number of effective interferometric pairs required to meet the temporal deformation inversion requirements of this monitoring sub-region. The conditions of a temporal baseline not exceeding 365 days and a vertical spatial baseline not exceeding 200 m are used for a preliminary evaluation of the suitability of the spatiotemporal baselines; when constructing SBAS-InSAR short-temporal-spatial baseline differential interferometric pair combinations for vegetation-covered areas, the constraints of a temporal baseline not exceeding 100 days and a spatial baseline not exceeding 150 m are further applied.

[0026] In one specific implementation, in step S2, the spatial overlay includes: selecting a reference SAR image from the SAR imagery, converting the land cover vector data to the coordinate system corresponding to the reference SAR imagery, reclassifying the land cover vector data, and overlaying the reclassification result with the landslide geological zone to obtain the boundary of the monitoring sub-zone. The reclassification is performed according to the influence of land cover type on radar scattering stability, classifying the original land cover type into at least one of bare land, forest land, grassland, and construction land; the landslide geological zone includes at least one of landslide backwall, landslide deposit, landslide impact zone, and relatively stable zone. By overlaying the reclassification result with the landslide geological zone, the spatial distribution relationship of different land cover types in the landslide geological unit is determined, and the boundary of the monitoring sub-zone is formed accordingly.

[0027] In another specific embodiment, the monitoring sub-area includes exposed bedrock area, landslide accumulation area, vegetation cover area and artificial building area, and is delineated according to the superposition relationship between landslide geological attributes and land cover type.

[0028] In this embodiment, the exposed bedrock area is defined by the back wall of the landslide, the bedrock exposure area, and the area with bare land cover; the landslide deposit area is defined by the distribution area of ​​the main body of the landslide and loose deposits as determined by the landslide geological survey data; the vegetation cover area is defined by the forest, grassland, shrub cover area, and vegetation cover area located within the landslide impact area; and the artificial building area is defined by the distribution area of ​​construction land, roads, villages, and other artificial structures.

[0029] Specifically, when the same spatial unit meets multiple delineation conditions simultaneously, it is first determined whether it belongs to a landslide accumulation zone based on the landslide's geological attributes. If it does not belong to a landslide accumulation zone, it is then classified into an exposed bedrock zone, a vegetation-covered zone, or an artificial building zone based on the land cover type. Through the above delineation rules, the monitoring sub-zones reflect both differences in land cover and differences in landslide geological attributes, providing a spatial basis for subsequent matching InSAR processing strategies.

[0030] In another specific embodiment, in step S3, the scatterer distribution characteristics include permanent scatterer candidate point distribution characteristics and distributed scatterer candidate point distribution characteristics. The permanent scatterer candidate point distribution characteristics are determined based on the pixel distribution that satisfies preset coherence conditions and preset intensity conditions. The distributed scatterer candidate point distribution characteristics are determined based on the statistical consistency between candidate pixels and neighboring pixels. The matching of the InSAR processing strategy includes: matching the PS-InSAR processing strategy to the exposed bedrock area and the artificial building area; matching the PS-InSAR combined with DS-InSAR processing strategy to the landslide deposit area; and matching the SBAS-InSAR processing strategy to the vegetation cover area.

[0031] Furthermore, the processing adaptation indicators include land cover type, average coherence coefficient, average intensity value, distribution characteristics of permanent scatterer candidate points, distribution characteristics of distributed scatterer candidate points, radar visibility index, and spatiotemporal baseline suitability. The distribution characteristics of permanent scatterer candidate points characterize the distribution of point-like stable scattering targets within the monitoring sub-area; the distribution characteristics of distributed scatterer candidate points characterize the distribution of area-like or regional scattering targets with consistent neighborhood statistical characteristics within the monitoring sub-area; the radar visibility index characterizes the degree to which the monitoring sub-area is affected by overlay, shadowing, perspective contraction, and low observation sensitivity, and is determined based on the area proportion and spatial distribution of low-confidence areas in the radar visibility mask within the monitoring sub-area; the spatiotemporal baseline suitability characterizes whether the SAR image assemblage meets the requirements for short spatiotemporal baseline interferometry processing. The judgment criteria for the processing adaptation indicators are determined based on the statistical distribution of quality indicators of each monitoring sub-area during the monitoring period, the quality indicators of the reference stable sample area, and the minimum effective point distribution requirements for InSAR processing.

[0032] Specifically, for monitoring sub-regions where the average coherence coefficient and average intensity value meet the requirements for permanent scatterer point selection, and the distribution of permanent scatterer candidate points meets the requirements for deformation calculation network construction, a PS-InSAR processing strategy is matched; for monitoring sub-regions where there are insufficient permanent scatterer candidate points but distributed scatterer candidate points meet the requirements for phase optimization and joint deformation calculation, a PS-InSAR combined with DS-InSAR processing strategy is matched; for monitoring sub-regions where the distribution of both permanent and distributed scatterer candidate points is insufficient, but the combination of short-temporal baseline interferometry pairs meets the requirements for temporal deformation inversion, a SBAS-InSAR processing strategy is matched. The above matching process is completed based on processing adaptation indicators, rather than by manually specifying the processing technology.

[0033] In another specific embodiment, when performing time-series deformation inversion on the exposed bedrock area and the artificial building area, permanent scatterer points are selected based on the average coherence coefficient in the coherence quality index and the average intensity value in the intensity quality index, and a PS-InSAR deformation solution network is constructed based on the permanent scatterer points.

[0034] Furthermore, the permanent scattering point is a point-like scattering target with stable phase characteristics and stable backscattering intensity during the monitoring period. The ground structures in exposed bedrock areas and artificial building areas are stable, and the radar backscattering characteristics change little over time, making them suitable for using permanent scattering points to construct spatial connections and perform temporal deformation calculations.

[0035] Specifically, in the PS-InSAR deformation calculation process, the phase sequence of the permanent scatterer points is associated with a deformation reference datum established by a unified reference point or a set of reference points. The temporal deformation and deformation rate of the permanent scatterer points are estimated, and the errors caused by atmospheric phase, orbital residuals, and topographic residuals are corrected or reduced to obtain the temporal deformation results of the exposed bedrock area and the artificial building area. In this invention, the temporal deformation is preferably the radar line-of-sight temporal deformation. When it is necessary to obtain the slope aspect deformation or three-dimensional deformation, the radar line-of-sight temporal deformation is projected and transformed or jointly calculated based on the ascending-descending orbit SAR image, slope aspect parameters, or external GNSS monitoring data.

[0036] In another specific embodiment, when performing time-series deformation inversion on the landslide deposit area, permanent scatterer points are first screened within the landslide deposit area, and a deformation calculation benchmark is established based on the permanent scatterer points within the landslide deposit area. Then, distributed scatterer points are screened according to the distribution characteristics of homogeneous scatterers, and joint deformation calculation is performed on the permanent scatterer points and the distributed scatterer points. The distributed scatterer points are determined by comparing the intensity statistical distribution, coherence coefficient change, and phase residual dispersion of candidate pixels with those of neighboring pixels.

[0037] Furthermore, the homogeneous scatterer distribution characteristics include the intensity statistical consistency, coherence stability, and time-series phase stability between candidate pixels and their neighboring pixels. Distributed scatterer points selected based on these homogeneous scatterer distribution characteristics are used to supplement sparse monitoring areas caused by insufficient permanent scatterer points within the landslide accumulation zone.

[0038] Furthermore, using the permanent scatterer points or the unified reference point or reference point set as a benchmark, phase optimization and deformation estimation are performed on the distributed scatterer points. The phase optimization is used to improve the phase quality of the distributed scatterer points, and the deformation estimation is used to obtain the deformation of the distributed scatterer points at different monitoring times. The deformation results of the permanent scatterer points and the distributed scatterer points are jointly constrained to obtain the time-series deformation results of the landslide accumulation area.

[0039] In one specific implementation, when performing temporal deformation inversion on the vegetation cover area, a short-spatial-temporal baseline differential interferometric pair combination is constructed based on the spatiotemporal baseline constraints, and SBAS-InSAR temporal deformation inversion is performed based on the short-spatial-temporal baseline differential interferometric pair combination.

[0040] Specifically, the spatiotemporal baseline constraints include temporal baseline constraints and spatial baseline constraints, which are determined based on the spatiotemporal baseline distribution of SAR images and the required number of effective interferometric pairs. Based on these temporal and spatial baseline constraints, short-spatial baseline differential interferometric pair combinations are constructed from SAR images within the monitoring period to reduce the impact of temporal and spatial incoherence in vegetation cover areas on deformation inversion.

[0041] Furthermore, the phase unwrapping of the short-spatial-time baseline differential interferometry pair combination is performed to establish the temporal deformation observation relationship. Combined with the deformation reference benchmark established by the unified reference point or reference point set, the atmospheric phase, orbital error and topographic residual are corrected or weakened to obtain the temporal deformation and deformation rate of the vegetation cover area.

[0042] In one specific implementation, in step S4, the unified reference point or reference point set is selected from the overlapping buffer zone of adjacent monitoring sub-regions, and the unified reference point or reference point set satisfies preset coherence conditions, preset intensity conditions, and preset temporal deformation stability conditions; the temporal deformation stability means that the absolute value of the deformation rate, the range of deformation fluctuation, or the phase residual of the candidate reference point during the monitoring period satisfies the preset stability conditions; the unified reference point or reference point set is a stable point or set of stable points selected from the overlapping buffer zone of adjacent monitoring sub-regions, used to establish a deformation reference benchmark commonly adopted by each monitoring sub-region.

[0043] In this embodiment, the buffer zone is a spatial region formed by extending a predetermined distance inward or outward along the boundary of the monitoring sub-region; the overlapping region of the buffer zones is the spatial overlap of the buffer zones of adjacent monitoring sub-regions. When three or more monitoring sub-regions are adjacent, the common overlapping region of the buffer zones of multiple monitoring sub-regions is used as a reference candidate region.

[0044] Specifically, within the reference candidate region, points with stable average coherence coefficients, stable average intensity values, and stable deformation changes during the monitoring period are selected as unified reference points. When a single stable point is insufficient to constrain all monitoring sub-regions, multiple stable points are selected to form a reference point set, and a unified deformation reference benchmark is established based on this reference point set. The time-series deformation inversion results of each monitoring sub-region are constrained to the unified deformation reference benchmark to reduce benchmark differences caused by independent inversions of different regions.

[0045] In one specific implementation, in step S5, the seamless spatial fusion includes: calculating the partitioned deformation residuals based on corresponding points, adjacent points, or grid cells within the overlapping areas of adjacent monitoring sub-regions; performing systematic error correction and regional network adjustment based on the partitioned deformation residuals; and sampling the corrected temporal deformation results onto a unified spatial resolution grid. The regional network adjustment uses the deviation correction amount of each monitoring sub-region as the unknown quantity, the deformation residuals of corresponding points, adjacent points, or grid cells within the overlapping areas as the observed quantity, and determines the observation weights based on the coherence quality index, intensity quality index, and inversion residuals. The seamless spatial fusion refers to sampling the deformation results of each monitoring sub-region after completing systematic error correction and regional network adjustment onto a unified spatial resolution grid, and determining the fusion result at the partition boundaries based on the quality index.

[0046] In this embodiment, the corresponding point is a monitoring point whose spatial position is consistent in the temporal deformation results of different monitoring sub-regions or falls into the same grid cell after spatial registration; the adjacent point is a monitoring point located near the boundary of adjacent monitoring sub-regions and whose spatial distance satisfies the proximity relationship; the partition deformation residual is the deformation difference between adjacent monitoring sub-regions at the corresponding point, adjacent point or grid cell.

[0047] Specifically, based on the partitioned deformation residuals, the overall offset and spatial trend deviation between each monitoring sub-region are identified, and systematic error correction is performed on the temporal deformation results of each monitoring sub-region. Monitoring points with abnormal residuals or abnormal quality indicators are eliminated or downweighted. Then, using the deformation residuals in the overlapping area as redundant observations, a regional network adjustment relationship is established to solve for the deviation correction amount of each monitoring sub-region. After the correction is completed, the temporal deformation results of each monitoring sub-region are sampled to a unified spatial resolution grid, and the fusion result is determined at the partition boundary based on the coherence quality index, intensity quality index, inversion residuals, and distribution of effective monitoring points to obtain the temporal deformation monitoring results under a unified reference benchmark for the entire landslide area. The outlier elimination refers to excluding monitoring points whose deformation residuals, inversion residuals, or quality indicators do not meet the preset conditions from the adjustment and fusion constraints, or reducing their weights. The systematic error correction includes correcting the overall offset and spatial trend deviation between different monitoring sub-regions based on the partitioned deformation residuals.

[0048] In a specific implementation, taking the Baoxing River East Basin in Fengtongzhai Township, Baoxing County, Sichuan Province as an example, the automatic landslide zoning processing method based on InSAR technology of the present invention is illustrated. This embodiment is only used to explain the technical solution of the present invention and does not constitute a limitation on the scope of protection of the present invention; in this embodiment, the study area includes the landslide back wall, landslide deposits, slope vegetation cover area, riverside roads, and village and town building areas. Multiple SAR images, orbital ephemeris data, digital elevation model data, land cover vector data, and landslide geological survey data are acquired during the continuous monitoring period of the study area. The landslide geological survey data includes the landslide boundary, landslide back wall boundary, landslide deposit boundary, landslide impact zone boundary, and the relatively stable zone determined by the survey.

[0049] In this embodiment, the method includes the following steps: Step 1: Data Acquisition and Preprocessing Multiple SAR images covering the study area were selected, and orbital ephemeris data corresponding to each SAR image were acquired. Based on the temporal and spatial baseline distribution of each SAR image, a SAR image located in the middle of the time series and with a relatively centered spatial baseline was selected as the reference SAR image. Orbital correction was performed on each SAR image based on the orbital ephemeris data, and the remaining SAR images were registered to the reference SAR image. Then, terrain phase was simulated using digital elevation model data, and the terrain phase was removed during interferometry processing to obtain a differential interferogram. After filtering the differential interferogram, the average coherence coefficient and average intensity value were calculated pixel-by-pixel for the monitoring period. In this embodiment, the average coherence coefficient is used as a coherence quality index to evaluate the phase stability of a pixel during the monitoring period; the average intensity value is used as an intensity quality index to evaluate the radar backscattering intensity and its stability of the pixel. The aforementioned quality indicators are used for subsequent monitoring sub-region processing strategy matching, scatterer screening, and unified reference point selection. Meanwhile, based on digital elevation model data, the incident angle and orbital direction of SAR images, and the slope aspect and gradient of the study area, a radar visibility mask is generated to identify overlapping areas, shadow areas, perspective contraction areas, and low observation sensitivity areas, and these areas are identified as low confidence areas. Step 2: Land cover reclassification and monitoring sub-zone delineation The land cover vector data was converted to the reference SAR image coordinate system, so that the land cover vector data, SAR image pixels, and landslide geological survey data were under the same spatial reference. Subsequently, the land cover vector data was reclassified, and the original land cover types were classified into four categories: bare land, forest land, grassland, and construction land. The reclassified land cover types are then spatially overlaid with the landslide backwall, landslide deposits, landslide impact zone, and relatively stable zone from the landslide geological survey data to automatically delineate four types of monitoring sub-zones. Specifically, areas with obvious bedrock exposure and bare land cover in the landslide backwall are designated as bare bedrock zones; areas containing the main body of the landslide and loose deposits are designated as landslide deposit zones; areas covered by forest, grassland, and shrubs are designated as vegetation zones; and areas with roads, village buildings, slope protection projects, and other construction land are designated as artificial building zones. When the same spatial unit meets multiple delineation conditions, it is first determined whether it belongs to a landslide deposit zone based on the landslide geological attributes. If it does not belong to a landslide deposit zone, it is then classified into a bare bedrock zone, a vegetation zone, or an artificial building zone based on the land cover type. Through this processing, the boundaries of the monitoring sub-zones corresponding to the SAR image space are obtained. Step 3: Calculating Adaptation Indicators and Matching InSAR Processing Strategies For each monitoring sub-region, the average coherence coefficient, average intensity value, distribution of permanent scatterer candidate points, and distribution of distributed scatterer candidate points are statistically analyzed. Based on the radar visibility mask, the area ratio and spatial distribution of overlapping areas, shadow areas, perspective contraction areas, and low observation sensitivity areas within the monitoring sub-region are statistically analyzed to determine radar visibility indicators. The temporal and spatial baseline suitability is determined based on the temporal and spatial baselines of the SAR image and the number of differential interferometric pairs satisfying preset temporal and spatial baseline constraints. The land cover type, average coherence coefficient, average intensity value, distribution characteristics of permanent scatterer candidate points, distribution characteristics of distributed scatterer candidate points, radar visibility indicators, and temporal and spatial baseline suitability of the monitoring sub-region are used together as processing adaptation indicators. In this embodiment, if the average coherence coefficient and average intensity value are relatively stable within a monitoring sub-region, and the permanent scatterer candidate points can form a continuous spatial distribution, then the monitoring sub-region is matched with the PS-InSAR processing strategy. Therefore, PS-InSAR processing strategies are matched for exposed bedrock areas and artificial structure areas. If the number of permanent scatterer candidate points in a monitoring sub-area is insufficient, but homogeneous scatterers with similar intensity statistical characteristics and phase stability exist in the neighborhood, then the monitoring sub-area is matched with a PS-InSAR combined with DS-InSAR processing strategy. Similarly, PS-InSAR combined with DS-InSAR processing strategies are matched for landslide deposit areas. If a monitoring sub-area is affected by vegetation cover, resulting in insufficient distribution of permanent and distributed scatterer points, but meets the conditions for constructing short-temporal baseline differential interferometry pairs, then the monitoring sub-area is matched with a SBAS-InSAR processing strategy. Therefore, the SBAS-InSAR processing strategy is matched to the vegetation-covered area. The matching process is completed automatically based on the processing adaptation index, rather than manually specifying different InSAR processing technologies. In this embodiment, for exposed bedrock areas and artificial building areas, pixels with an average coherence coefficient of not less than 0.40 and an average intensity value not less than the 75th percentile of the intensity statistical distribution of the study area are selected as permanent scatterer candidate points. For landslide deposit areas, pixels with an average coherence coefficient of not less than 0.20 and whose intensity statistical distribution and phase residual satisfy the homogeneity condition in the 7×7 neighborhood are selected as distributed scatterer candidate points. For vegetation-covered areas, SAR image pairs with a temporal baseline of not more than 100 days and a spatial baseline of not more than 150m are constructed into short-temporal-spatial baseline differential interferometric pair combinations. The above thresholds are determined based on the statistical results of SAR image quality in the study area, the reference stable area samples, and the effective point distribution requirements for time-series deformation inversion. For low-confidence areas marked by radar visibility masks, the matching weight of the corresponding pixels is reduced during the processing strategy matching and scatterer screening process to reduce the impact of overlay, shadows, perspective shrinkage, and low observation sensitivity on the processing strategy matching results. Step 4: Selecting a unified reference point or set of reference points Buffer zones are established at the adjacent boundaries of exposed bedrock areas, landslide accumulation areas, vegetated areas, and artificial construction areas, and the overlapping areas of the buffer zones of adjacent monitoring sub-areas are obtained. For areas adjacent to multiple monitoring sub-areas, the common overlapping part of multiple buffer zones is taken as a reference candidate area. Within the reference candidate area, points with stable average coherence coefficients, stable average intensity values, and stable temporal deformation changes during the monitoring period are selected as candidate reference points. If a single candidate reference point can simultaneously constrain the deformation inversion of each monitoring sub-area, it is determined as a unified reference point; if a single candidate reference point cannot cover all monitoring sub-areas, multiple candidate reference points are selected to form a reference point set, and a unified deformation reference benchmark is established based on the weighted average deformation of the reference point set. In this embodiment, the unified reference point or reference point set is located in a relatively stable area outside the main deformation zone of the landslide and is spatially associated with at least two adjacent monitoring sub-areas. When performing time-series deformation inversion in each monitoring sub-region, the deformation reference benchmark established by the unified reference point or reference point set is used as a constraint. In this embodiment, a buffer zone is formed by extending 30m outwards and inwards from the boundary of each monitoring sub-region, and the overlapping area of ​​the buffer zones of adjacent monitoring sub-regions is used as the reference candidate area. Within the reference candidate area, points with an average coherence coefficient of not less than 0.45, an average intensity value of not less than the 75th percentile of the intensity statistical distribution of the study area, and an absolute value of deformation rate of not more than 2mm / a during the monitoring period are selected as unified reference points. When a single reference point cannot cover the entire monitoring sub-region, multiple points that meet the above conditions are selected to form a reference point set. Step 5: Temporal Deformation Inversion of Exposed Bedrock Areas and Artificial Building Areas For exposed bedrock areas and artificial building areas, temporal deformation inversion is performed according to the PS-InSAR processing strategy. Specifically, permanent scatterer points are screened using the average coherence coefficient and average intensity value. The screened permanent scatterer points include stable scattering points on exposed bedrock surfaces, stable structural points on road slopes, building corner points, and other strong scattering target points. Subsequently, a PS-InSAR deformation calculation network is constructed based on the permanent scatterer points, and the unified reference point or reference point set determined in S4 is used as the deformation reference benchmark. During the deformation calculation process, the phase sequence of the permanent scatterer points is modeled to reduce the influence of atmospheric phase, orbital residuals, and topographic residuals on the deformation results, thereby obtaining the temporal deformation and deformation rate at each monitoring time in the exposed bedrock area and artificial building area. Step Six: Temporal Deformation Inversion of the Landslide Accumulation Area For landslide accumulation areas, temporal deformation inversion is performed using a PS-InSAR combined with DS-InSAR processing strategy. Specifically, permanent scatterer points with satisfactory stability are first selected within the landslide accumulation area, and a local deformation solution network for the landslide accumulation area is established based on these permanent scatterer points. Then, distributed scatterer points are selected based on the homogeneous scatterer distribution characteristics, and joint deformation solutions are performed on the permanent scatterer points and the distributed scatterer points under the deformation reference constraints established by the unified reference point or reference point set. The homogeneous scatterer distribution characteristics include the intensity statistical consistency, coherence stability, and temporal series phase stability between candidate pixels and neighboring pixels. For the selected distributed scatterer points, phase optimization is performed using the permanent scatterer points or the unified reference point set as a reference to reduce the phase noise of the distributed scatterer points. After phase optimization, joint deformation calculation is performed on permanent scatterer points and distributed scatterer points to obtain the temporal deformation and deformation rate of the landslide deposit area. Through this process, the deformation monitoring results can be supplemented by distributed scatterer points in areas of the landslide deposit area that were originally difficult to cover by permanent scatterer points alone. Step 7: Temporal Deformation Inversion of Vegetation Cover Areas For vegetated areas, temporal deformation inversion is performed according to the SBAS-InSAR processing strategy. Specifically, based on the temporal and spatial baseline relationships between SAR images, SAR image pairs that meet short-term spatiotemporal baseline constraints are selected to construct short-term spatiotemporal baseline differential interferometric pair combinations. Subsequently, differential interferometric processing and phase unwrapping are performed on the short-term spatiotemporal baseline differential interferometric pair combinations to establish the temporal deformation observation relationship of the vegetated area. During the temporal deformation inversion process, the deformation results of the vegetated area are constrained by the deformation reference benchmark established by a unified reference point or reference point set, and atmospheric phase, orbital error, and topographic residuals are corrected or weakened to obtain the temporal deformation and deformation rate of the vegetated area at each monitoring time. Step 8: Partition Result Correction and Spatial Fusion The temporal deformation results of exposed bedrock areas, artificial structure areas, landslide deposit areas, and vegetation-covered areas are uniformly transformed to the same coordinate system and correspond to the deformation reference benchmark established by a unified reference point or reference point set. Within the overlapping buffer zones of adjacent monitoring sub-regions, deformation results of corresponding points, adjacent points, or grid cells are extracted, and the partition deformation residuals between different monitoring sub-regions are calculated. Corresponding points refer to monitoring points with consistent spatial locations in different partition processing results or falling into the same grid cell after spatial registration; adjacent points refer to monitoring points located near the boundaries of adjacent partitions and whose spatial distance satisfies the proximity relationship. Based on the partition deformation residuals, systematic error correction is performed on the temporal deformation results of each monitoring sub-region. Systematic error correction includes correcting the overall offset and spatial trend deviation between different partition results. For monitoring points with significantly abnormal deformation residuals or quality indicators, outlier removal or weight reduction is performed. Subsequently, using the partition deformation residuals within the overlapping area as redundant observations, a regional network adjustment relationship is established, and the deviation correction amount for each monitoring sub-region is solved. The deviation correction is applied to the temporal deformation results of the corresponding monitoring sub-region to maintain the continuity of deformation results in the boundary region of adjacent monitoring sub-regions. In this embodiment, deformation residuals are calculated for corresponding points, adjacent points, or grid cells in the overlapping area of ​​adjacent monitoring sub-regions. When the deformation residual exceeds three times the standard deviation of the residual in the overlapping area, the corresponding monitoring point is treated as an outlier and removed or downweighted. The regional network adjustment uses the deviation correction of each monitoring sub-region as the unknown quantity and the deformation residual in the overlapping area as the observed quantity, and determines the observation weights based on the average coherence coefficient, average intensity value, and inversion residual. After the adjustment is completed, the temporal deformation results of each monitoring sub-region are sampled to a 30m×30m uniform spatial resolution grid to obtain the temporal deformation monitoring results under the unified reference benchmark of the entire landslide area. For grid cells with multiple sets of deformation results at the boundary of the zone, the fusion result is determined based on the coherence quality index, intensity quality index, inversion residual, and distribution of effective monitoring points to form the temporal deformation monitoring results under the unified reference benchmark of the entire landslide area.

[0050] Specifically, after processing, the following results are output: Temporal deformation results of permanent scatterer points in exposed bedrock areas and artificial building areas; joint temporal deformation results of permanent and distributed scatterer points in landslide deposit areas; SBAS-InSAR temporal deformation results in vegetated areas; deformation rate map under a unified reference benchmark for the entire landslide area; temporal deformation map under a unified reference benchmark for the entire landslide area; deformation results of a unified grid after seamless spatial fusion of partition boundaries; Through the above specific embodiments, a complete processing flow can be realized from data acquisition, automatic partitioning, processing strategy matching, unified reference benchmark establishment, partition temporal deformation inversion to partition result fusion.

[0051] In one specific implementation, to verify the improvement effect of the method of the present invention on monitoring point coverage in complex landslide areas, deformation extraction accuracy in low-coherence areas, and continuity of zoning results, a typical experimental area in western Sichuan was selected as the study area. Eighty-seven Sentinel-1 ascending-orbit SAR images acquired between January 2018 and February 2021 were used for comparative verification. The SAR images were C-band VV polarimetric images with a relative orbit number of 99 and a strip frame number of 405. The control processing and the processing of the present invention used the same batch of SAR images, orbital ephemeris data, digital elevation model data, land cover vector data, and landslide geological survey data.

[0052] This embodiment sets up two sets of comparative verification. The first set of comparative verification is used to compare the differences between traditional unified PS-InSAR processing and the automatic partitioning processing of the present invention in terms of the number of monitoring points, coverage, monitoring capability in low coherence areas, deformation consistency between regions, and deformation inversion accuracy; the second set of comparative verification is used to compare the changes in partition boundary deformation residuals and boundary continuity before and after spatial fusion using the method of the present invention.

[0053] In the first comparative verification, the control group did not delineate monitoring sub-regions or perform processing strategy matching. The entire landslide area was uniformly processed using PS-InSAR, and pixels with an average coherence coefficient of not less than 0.40 and an average intensity value not lower than the 75th percentile of the intensity statistical distribution of the study area were selected as candidate permanent scatterers. The invention group processed the area according to steps S1 to S5, dividing the study area into exposed bedrock, landslide deposit, vegetation cover, and artificial structure areas, and matching PS-InSAR, PS-InSAR combined with DS-InSAR, and SBAS-InSAR processing strategies respectively. To objectively compare the spatial consistency of deformation results obtained from traditional unified PS-InSAR processing and the automatic partitioning processing of this invention across different scattering characteristic regions, the exposed bedrock, landslide deposit, vegetation cover, and artificial structure areas automatically delineated by the method of this invention were used as a unified partitioning benchmark, and the inter-regional deformation residuals of the two processing results were calculated respectively.

[0054] For the traditional unified PS-InSAR processing results, the deformation monitoring results are divided according to the boundaries of the four types of monitoring sub-regions mentioned above. In the overlapping area of ​​the 30m buffer zone at the boundary of adjacent monitoring sub-regions, the standard deviation of the deformation difference between the same point or adjacent points is calculated as the inter-regional deformation residual of the traditional unified PS-InSAR processing results.

[0055] For the automatic partitioning processing results of the present invention, after completing the system error correction, regional network adjustment and spatial seamless fusion in step S5, the standard deviation of the deformation difference between the same monitoring sub-region boundary and the same 30m buffer overlap area is statistically analyzed as the inter-regional deformation residual of the automatic partitioning processing results of the present invention.

[0056] Although the traditional unified PS-InSAR processing does not explicitly delineate monitoring sub-regions, the same four types of monitoring sub-region boundaries as the method of this invention are used as a unified evaluation benchmark in the experimental evaluation stage to ensure that the deformation residuals between regions of the two processing results have the same calculation area and calculation caliber.

[0057] Figure 2 , Figure 3 The spatial distribution of monitoring points for the two sets of processing results is shown. Figure 2 It can be seen that the control group had relatively concentrated monitoring points in the exposed bedrock area and the artificial structure area, but sparse monitoring points and monitoring gaps existed in the landslide deposit area and the vegetation cover area. Figure 3 It can be seen that the present invention group has improved the number of monitoring points and the continuity of spatial distribution in the above-mentioned areas by introducing distributed scatterer points for landslide accumulation areas and constructing short-temporal baseline differential interferometry pairs for vegetation-covered areas.

[0058] like Figure 4 As shown, the control group obtained 7,000, 1,100, 700 and 2,800 effective monitoring points in the four types of monitoring sub-areas, respectively, for a total of 11,600 effective monitoring points in the entire landslide area; the invention group obtained 7,000, 5,200, 4,100 and 2,800 effective monitoring points in the four types of monitoring sub-areas, respectively, for a total of 19,100 effective monitoring points in the entire landslide area, an increase of approximately 65% ​​compared to the control group.

[0059] This embodiment uses the entire landslide area boundary as the statistical scope, dividing the entire landslide area into a uniform spatial resolution grid of 30m × 30m, with each grid corresponding to an area of ​​0.0009 km². Grids that yield effective temporal deformation results under the corresponding InSAR processing strategy for each monitoring sub-region are defined as valid grids, provided they contain at least one valid monitoring point that meets the corresponding scatterer screening criteria. The effective area coverage rate is determined based on the proportion of valid grids within the corresponding area to the total number of grids in that area.

[0060] In one implementation, the entire landslide area comprises 10,234 grid cells, covering an area of ​​9.21 km². Specifically, the exposed bedrock area comprises 4,767 grid cells, covering an area of ​​4.29 km², accounting for 46.58% of the total landslide area; the landslide deposit area comprises 1,733 grid cells, covering an area of ​​1.56 km², accounting for 16.93% of the total landslide area; the vegetation-covered area comprises 967 grid cells, covering an area of ​​0.87 km², accounting for 9.45% of the total landslide area; and the artificial structure area comprises 2,767 grid cells, covering an area of ​​2.49 km², accounting for 27.04% of the total landslide area.

[0061] Statistical analysis showed that the control group had 4,528 effective grid cells in the exposed bedrock area, landslide deposit area, vegetated area, and artificial building area, respectively, with corresponding effective area coverage rates of 94.99%, 35.03%, 12.00%, and 98.01%. Figure 5 The percentages are rounded down to 95%, 35%, 12%, and 98%. The control group had 7,963 effective grids throughout the landslide area, with an effective area coverage rate of 77.81%, approximately 78%.

[0062] The number of effective grid cells in the exposed bedrock area, landslide accumulation area, vegetation cover area, and artificial construction area of ​​this invention are 4528, 1473, 696, and 2712, respectively, with corresponding effective area coverage rates of 94.99%, 85.00%, 71.98%, and 98.01%. Figure 5 The percentages are rounded down to 95%, 85%, 72%, and 98%. The invention group has a total of 9409 effective grids throughout the entire landslide area, with an effective area coverage rate of 91.94%, approximately 92%, which is 14.13 percentage points higher than the control group.

[0063] The effective area coverage rate of the entire landslide area was calculated by summing the number of effective grid cells in the four types of monitoring sub-regions and the total number of grid cells in the entire landslide area. Since each grid cell has the same spatial resolution and corresponding area, the effective area coverage rate of the entire landslide area can be obtained by summing the number of effective grid cells in the four types of monitoring sub-regions.

[0064] Figure 6 The deformation monitoring results in low-coherence areas are represented by a dual ordinate system, where the left ordinate represents the number of effective monitoring points in the low-coherence area (in units), and the right ordinate represents the percentage of void area in the low-coherence area. The low-coherence area includes landslide deposits and vegetation cover. Traditional unified PS-InSAR processing yielded 1800 effective monitoring points in the low-coherence area, with a void area percentage of 35%. The automatic ordinate system of this invention yielded 9300 effective monitoring points in the low-coherence area, an increase of approximately 417% compared to traditional unified PS-InSAR processing, while reducing the void area percentage to 8%, a decrease of 27 percentage points, or approximately 77%. The void area refers to a continuous area without effective monitoring results identified through grid connectivity analysis; scattered invalid grids are not included in the void area calculation. The void area percentage is the ratio of the sum of the areas of all continuous void areas within the low-coherence area to the total area of ​​the low-coherence area. Therefore, the percentage of void area and the effective area coverage rate are calculated using different statistical methods and are not directly complementary.

[0065] Based on the above unified zoning benchmark, the inter-regional deformation residuals are calculated as follows: Figure 7 As shown, the inter-regional deformation residual of the traditional unified PS-InSAR processing result is 9.6 mm, while the inter-regional deformation residual of the automatic partitioning processing result of this invention is 2.1 mm, a reduction of approximately 78% compared to the traditional unified PS-InSAR processing. This result indicates that although the traditional unified PS-InSAR processing does not explicitly delineate monitoring sub-regions, there are still significant differences in deformation results between regions with different scattering characteristics. This invention, through automatic delineation of monitoring sub-regions, matching of processing strategies, unified deformation reference benchmarks, and seamless spatial fusion, reduces the systematic differences between regions with different scattering characteristics and improves the spatial consistency of deformation results across the entire landslide area. Using the cumulative deformation along the radar line of sight at six global navigation satellite system monitoring points in the experimental area as a reference, the root mean square error (RMSE) between the two sets of processed results and the field monitoring results is calculated. Figure 7 As shown, the RMSE of the control group was 7.1 mm, while the RMSE of the present invention group was 3.4 mm, which was about 52% lower than that of the control group.

[0066] Table 1. Comprehensive Comparison of Traditional Unified PS-InSAR and the Automatic Partitioning Processing of the Invention

[0067] Table 1 shows that the automatic zoning processing of this invention, while maintaining monitoring capabilities in exposed bedrock and artificial construction areas, increases the number and coverage of effective monitoring points in landslide accumulation and vegetation-covered areas. Using the same monitoring sub-region boundaries and buffer zones as a unified evaluation benchmark, the automatic zoning processing of this invention reduces the inter-regional deformation residual from 9.6 mm to 2.1 mm and the GNSS verification point RMSE from 7.1 mm to 3.4 mm. These results indicate that matching the corresponding InSAR processing strategy to the processing adaptation index of the monitoring sub-region, and further implementing unified reference benchmark constraints and seamless spatial fusion, can reduce the loss of deformation information in low-coherence areas when using single PS-InSAR processing, and improve the spatial consistency of deformation results between areas with different scattering characteristics.

[0068] The second set of comparative verification uses the independent temporal deformation results of each monitoring sub-region obtained in step S4 as the result before fusion, and the result after completing systematic error correction, regional network adjustment, and unified spatial resolution grid sampling in step S5 as the result after fusion. In this embodiment, the standard deviation of the deformation difference between corresponding points or adjacent points in the overlapping area of ​​adjacent monitoring sub-regions is used as the partition boundary deformation residual. The partition boundary deformation residual is used to evaluate the systematic differences of the independent deformation inversion results of each monitoring sub-region obtained in step S4 of the method of the present invention at the partition boundary, and the correction effect of step S5 on the systematic differences. Since the InSAR processing strategies used by each monitoring sub-region are different, even if each temporal deformation result has been constrained to a unified deformation reference benchmark, adjacent monitoring sub-regions before fusion still have partition boundary deformation residuals in the buffer overlapping area.

[0069] Figure 8 Before fusion, the residual deformation of the partition boundary, the root mean square error of deformation in the overlapping area, and the standard deviation of deformation in the overlapping area were 8.4 mm, 5.2 mm, and 4.7 mm, respectively. After fusion, they decreased to 2.3 mm, 2.0 mm, and 2.1 mm, respectively. Figure 9 In the process, the number of boundary jump pixels decreased from 268 before fusion to 49 after fusion.

[0070] Table 2: Overall Comparison of Effects Before and After Partition Merging

[0071] Table 2 shows that step S5, through systematic error correction and regional network adjustment based on the zonal deformation residuals, reduced the zonal boundary deformation residuals of adjacent monitoring sub-regions from 8.4 mm to 2.3 mm, the RMSE of overlapping area deformation from 5.2 mm to 2.0 mm, and the standard deviation of overlapping area deformation from 4.7 mm to 2.1 mm, while also reducing the number of boundary jump pixels from 268 to 49. These results demonstrate that step S5 can correct the systematic deviations between the independent inversion results of each monitoring sub-region, improving the continuity and consistency of the time-series deformation monitoring results across the entire landslide area at the zonal boundaries.

[0072] Based on the above two sets of comparative verifications, the method of the present invention, through the combination of monitoring sub-region delineation, processing strategy matching, unified deformation reference benchmark and seamless spatial integration, has achieved corresponding improvements in the number of effective monitoring points, monitoring coverage, monitoring capability in low coherence areas, degree of conformity of on-site monitoring results and continuity of partition boundaries, thereby providing data support for the technical solution and technical effects of the present invention.

[0073] In summary, this embodiment has at least the following technical effects: By acquiring SAR imagery, orbital ephemeris data, digital elevation model data, land cover vector data, and landslide geological survey data for the monitoring period in the study area, and by preprocessing and calculating quality indicators for the SAR imagery, a unified data foundation is provided for subsequent zoning, scatterer screening, and temporal deformation inversion. Compared to methods that rely solely on SAR imagery for unified processing, this invention combines land cover information, landslide geological information, and InSAR quality indicators, improving the completeness and relevance of preliminary data processing for deformation monitoring in landslide areas. This invention spatially overlays land cover vector data with landslide geological zoning data from landslide geological surveys, automatically delineating exposed bedrock areas, landslide deposit areas, vegetation cover areas, and artificial structure areas. This allows the monitoring sub-regions to reflect both differences in surface cover types and differences in the internal geological structure of the landslide body, avoiding the inaccurate zoning problems caused by dividing areas solely by land cover type or solely by landslide geological boundaries. This provides a more reasonable spatial unit for subsequent InSAR processing strategy matching. This invention forms processing adaptation indicators based on the land cover type, coherence quality index, intensity quality index, scatterer distribution characteristics, radar visibility index, and spatiotemporal baseline suitability of each monitoring sub-area. Based on the processing adaptation indicators, it matches the corresponding InSAR processing strategy and scatterer screening rules, avoiding the subjectivity and instability caused by directly selecting processing technology based on human experience, and improving the automation level and reusability of the landslide InSAR processing process. This invention employs a PS-InSAR processing strategy for exposed bedrock areas and artificial building areas, a PS-InSAR combined with DS-InSAR processing strategy for landslide deposit areas, and a SBAS-InSAR processing strategy for vegetated areas. This allows different scattering characteristic areas to adopt time-series deformation inversion methods that are adapted to their ground features and quality indicators, thereby improving the effective monitoring point density in complex landslide areas and reducing the problem of missing deformation information in low-coherence areas and loose deposit areas. By selecting a unified reference point or set of reference points in the overlapping area of ​​the buffer zones of adjacent monitoring sub-regions, and establishing a deformation reference benchmark based on the unified reference point or set of reference points, the time-series deformation inversion results of each monitoring sub-region are under a unified reference benchmark, which reduces the risk of benchmark offset, inconsistent results and abrupt boundary deformation after independent calculation of different partitions. This invention performs spatial registration, systematic error correction, outlier removal, and regional network adjustment on the temporal deformation results of each monitoring sub-region, and samples the corrected results to a unified spatial resolution grid. This can reduce the systematic bias between different processing strategies, improve the continuity and consistency of deformation results at the boundary of the partition, and finally obtain the temporal deformation monitoring results under a unified reference benchmark for the entire landslide area. This invention forms a closed-loop process from data acquisition and preprocessing, automatic delineation of monitoring sub-regions, automatic matching of processing strategies, construction of a unified reference benchmark, partitioned temporal deformation inversion to seamless spatial fusion, which can improve the automation level, monitoring integrity and result reliability of InSAR deformation monitoring in complex landslide areas.

[0074] It is understood that the above embodiments are merely exemplary implementations used to illustrate the principles of the present invention, and the present invention is not limited thereto. For those skilled in the art, various modifications and improvements can be made without departing from the spirit and essence of the present invention, and these modifications and improvements are also considered to be within the scope of protection of the present invention.

Claims

1. An automatic landslide zoning processing method based on InSAR technology, characterized in that, Includes the following steps: S1: Acquire SAR images, orbital ephemeris data, digital elevation model data, land cover vector data, and landslide geological survey data for the monitoring period in the study area; S2: Based on the orbital ephemeris data and the digital elevation model data, the SAR image is preprocessed, and the coherence quality index and intensity quality index for zoning processing are calculated; the land cover vector data is spatially overlaid with the landslide geological zoning in the landslide geological survey data to automatically delineate monitoring sub-zones; a radar visibility mask is generated based on the observation geometry of the digital elevation model data and the SAR image, and the radar visibility index is determined according to the distribution of the radar visibility mask in each monitoring sub-zone; and the spatiotemporal baseline suitability is determined according to the temporal and spatial baseline distribution of the SAR image. S3: Based on the land cover type, coherence quality index, intensity quality index, scatterer distribution characteristics, radar visibility index and spatiotemporal baseline suitability of each monitoring sub-region, a processing adaptation index is formed for each monitoring sub-region, and the corresponding InSAR processing strategy and scatterer screening rules are matched based on the processing adaptation index. S4: Based on the InSAR processing strategy and scatterer screening rules obtained by matching each monitoring sub-region, and combined with the deformation reference benchmark established by the unified reference point or reference point set, perform time-series deformation inversion on each monitoring sub-region to obtain the time-series deformation results of each monitoring sub-region. S5: Spatial registration, system error correction, outlier removal and spatial seamless fusion are performed on the temporal deformation results of each monitoring sub-region to obtain the temporal deformation monitoring results under a unified reference benchmark for the entire landslide area. The temporal deformation monitoring results include the radar line-of-sight temporal deformation.

2. The automatic landslide zoning processing method based on InSAR technology according to claim 1, characterized in that, In step S2, the preprocessing includes track correction, image registration, differential interferometry, filtering, and quality index calculation. The quality index calculation includes calculating the average coherence coefficient and average intensity value of the pixels.

3. The automatic landslide zoning processing method based on InSAR technology according to claim 1, characterized in that, In step S2, the spatial overlay includes: selecting a reference SAR image from the SAR image, converting the land cover vector data to the coordinate system corresponding to the reference SAR image, reclassifying the land cover vector data, and overlaying the reclassification result with the landslide geological zone to obtain the boundary of the monitoring sub-zone.

4. The automatic landslide zoning processing method based on InSAR technology according to claim 3, characterized in that, The monitoring sub-areas include exposed bedrock areas, landslide accumulation areas, vegetation-covered areas, and artificial structure areas, and are delineated according to the superposition relationship between landslide geological attributes and land cover types.

5. The automatic landslide zoning processing method based on InSAR technology according to claim 4, characterized in that, In step S3, the scatterer distribution features include permanent scatterer candidate point distribution features and distributed scatterer candidate point distribution features. The permanent scatterer candidate point distribution features are determined based on the pixel distribution that satisfies preset coherence and intensity conditions. The distributed scatterer candidate point distribution features are determined based on the statistical consistency between candidate pixels and neighboring pixels. The matching of the InSAR processing strategy includes: when the permanent scatterer candidate point distribution in the exposed bedrock area or the artificial building area satisfies the deformation solution network construction conditions, the PS-InSAR processing strategy is matched; when the permanent scatterer candidate points and distributed scatterer candidate points in the landslide deposit area satisfy the joint solution conditions, the PS-InSAR combined with DS-InSAR processing strategy is matched; when the vegetation cover area satisfies the short-spatial-temporal baseline differential interferometry pair construction conditions, the SBAS-InSAR processing strategy is matched.

6. The automatic landslide zoning processing method based on InSAR technology according to claim 5, characterized in that, The coherence quality index includes the average coherence coefficient, and the intensity quality index includes the average intensity value. When performing time-series deformation inversion on the exposed bedrock area and the artificial building area, permanent scatterer points are selected based on the average coherence coefficient and the average intensity value, and a PS-InSAR deformation solution network is constructed based on the permanent scatterer points.

7. The automatic landslide zoning processing method based on InSAR technology according to claim 5, characterized in that, When performing time-series deformation inversion on the landslide deposit area, permanent scatterer points are first selected within the landslide deposit area, and a local deformation solution network is established based on the permanent scatterer points. Then, distributed scatterer points are selected according to the distribution characteristics of homogeneous scatterers. Under the deformation reference datum constraint established by the unified reference point or reference point set, joint deformation solution is performed on the permanent scatterer points and the distributed scatterer points.

8. The automatic landslide zoning processing method based on InSAR technology according to claim 5, characterized in that, When performing temporal deformation inversion on the vegetation cover area, a short-spatial-temporal baseline differential interferometric pair combination is constructed based on the spatiotemporal baseline constraints, and SBAS-InSAR temporal deformation inversion is performed based on the short-spatial-temporal baseline differential interferometric pair combination.

9. The automatic landslide zoning processing method based on InSAR technology according to claim 1, characterized in that, In step S4, the unified reference point or reference point set is selected from the overlapping area of ​​the buffer zone of adjacent monitoring sub-regions, and the unified reference point or reference point set satisfies the preset coherence condition, preset intensity condition and preset temporal deformation stability condition.

10. The automatic landslide zoning processing method based on InSAR technology according to any one of claims 1 to 9, characterized in that, In step S5, the spatial seamless fusion includes: calculating the partition deformation residual based on the corresponding points, adjacent points or grid cells in the buffer overlap area of ​​adjacent monitoring sub-regions, performing system error correction and regional network adjustment based on the partition deformation residual, and sampling the corrected temporal deformation result to a unified spatial resolution grid.

Images