A landslide surface deformation spatiotemporal characteristic interpolation method, a storage medium and an equipment

Abstract

This invention discloses a method, storage medium, and device for interpolating spatiotemporal features of landslide surface deformation, belonging to the field of landslide disaster prediction and forecasting technology. The method includes the following steps: acquiring raw InSAR point cloud data using SBAS-InSAR technology; cleaning the data using the DBSCAN clustering algorithm that integrates spatial and deformation features; adaptively weighting and correcting the InSAR points using high-precision GNSS data based on KDTree and Gaussian kernel functions to eliminate system bias; constructing a heterogeneous fully connected graph containing geometric coordinates and environmental factors; and constructing a channel attention-enhanced physical guided graph transformation network SE-PG-GT, introducing a slope-constrained physical consistency loss function, and training it under a semi-supervised transductive learning framework to reconstruct the deformation field from sparse discrete points to a continuous grid. This invention effectively integrates multi-source data, providing reliable data support for accurate landslide early warning.

Description

Technical Field

[0001] This invention relates to the field of landslide disaster prediction and forecasting technology, and in particular to a method, storage medium and device for interpolating the spatiotemporal characteristics of landslide surface deformation. Background Technology

[0002] Landslides, a globally widespread geological hazard, often cause severe casualties and property damage due to their suddenness and destructiveness. Establishing an efficient and accurate landslide monitoring and early warning system is crucial for identifying risks in advance and guiding disaster prevention and mitigation. Surface deformation is a direct indicator of landslide stability evolution. Obtaining a high-precision, continuous planar deformation field covering the entire landslide area can accurately locate active deformation zones, providing core data for landslide control engineering deployment and disaster early warning.

[0003] Currently, methods for obtaining landslide surface displacement fields mainly rely on spatial interpolation of discrete monitoring points, commonly using traditional geostatistical models such as inverse distance weighted (IDW) and Kriging. However, landslide deformation is controlled by multiple factors, including topography, geological structure, and hydrological conditions, exhibiting significant spatial non-uniformity and anisotropy. Traditional interpolation methods are typically based on the isotropic stationarity assumption, making it difficult to effectively characterize nonlinear deformation gradients in complex mountain environments. This often results in overly smooth reconstructions, failing to preserve local deformation details and boundary features. Furthermore, relying solely on expensive Global Navigation Satellite System (GNSS) monitoring stations for interpolation is limited by deployment costs and the sparsity of monitoring points, making it difficult to accurately reflect the overall spatial motion morphology of landslides. Synthetic Aperture Radar Interferometry (InSAR) technology, with its advantages of all-weather operation, wide coverage, and high spatial resolution, has become an important means of monitoring micro-deformation of the landslide surface. However, InSAR technology is affected by atmospheric delay, spatiotemporal decoherence, and phase unwrapping errors. Its observation results often show discrete point cloud distributions, and it is prone to data holes and long-wavelength system bias in areas with dense vegetation or drastic terrain.

[0004] In summary, the existing landslide surface displacement reconstruction technology has the following main defects: (1) The traditional interpolation model is based on the isotropic assumption, which ignores the anisotropic and nonlinear characteristics of landslide deformation, making it difficult to accurately reconstruct the deformation gradient and details under complex terrain; (2) The sparse GNSS points cannot reflect the overall shape, and InSAR is prone to systematic bias and has voids. Existing methods lack high-precision physical fusion of heterogeneous data, making it difficult to balance spatial continuity and monitoring accuracy. Summary of the Invention

[0005] The purpose of this invention is to propose a method, storage medium, and device for interpolating the spatiotemporal characteristics of landslide surface deformation, thereby solving the technical problem that existing methods lack high-precision physical fusion of heterogeneous data and are difficult to balance spatial continuity and monitoring accuracy.

[0006] Specifically, the present invention provides a method for interpolating the spatiotemporal characteristics of landslide surface deformation, comprising the following steps: S1. Acquire SAR images of the target landslide body, process them using SBAS-InSAR technology, extract the displacement time series of the landslide body surface, and obtain the original InSAR monitoring point cloud data. S2. The original InSAR monitoring point cloud data is cleaned using the DBSCAN clustering algorithm that integrates spatial and deformation features to identify and remove noise points and outliers, resulting in a cleaned InSAR dataset. S3. Obtain observation data from GNSS stations in the study area, and use KDTree to construct local neighborhood indexes for each GNSS station. Based on the Gaussian kernel function, perform adaptive weighted correction on InSAR points falling into the neighborhood to eliminate systematic bias in InSAR data, realize spatial fusion of multi-source data, and obtain the corrected InSAR point dataset. S4. Construct a heterogeneous fully connected graph containing geometric coordinates and multi-source environmental factors, and use the corrected InSAR point data and the grid points to be interpolated as graph nodes to form graph structure data. S5. Construct and train the Physical Guided Graph Transformation Network (SE-PG-GT) model based on channel attention enhancement. Input the graph structure data into the SE-PG-GT model and learn the mapping relationship from sparse discrete points to continuous grids under the semi-supervised transductive learning framework. S6. Use the trained SE-PG-GT model to predict all nodes in the map and output a high-precision, physically constrained, continuous planar deformation field covering the entire landslide body.

[0007] A storage medium storing instructions and data for implementing a landslide surface deformation spatiotemporal feature interpolation method.

[0008] A landslide surface deformation spatiotemporal feature interpolation device includes: a processor and the storage medium; the processor loads and executes instructions and data in the storage medium to implement a landslide surface deformation spatiotemporal feature interpolation method.

[0009] The beneficial effects provided by this invention are: it achieves a deep integration of physical mechanisms and data-driven approaches, effectively improving the geological rationality of the interpolation results. Unlike traditional IDW or Kriging interpolation methods, which are usually based on isotropic spatial assumptions, the SE-PG-GT model proposed in this invention introduces a physical consistency loss function and uses prior knowledge of terrain slope to apply soft constraints to deformation prediction in areas without data. This effectively suppresses spurious deformation noise in flat areas and ensures that the generated deformation field strictly follows the physical laws of landslides driven by gravitational potential energy.

[0010] By using the channel attention mechanism (SE-Block) to dynamically filter key disaster-causing factors such as slope and aspect, the model can accurately capture nonlinear deformation gradients and local details under complex terrain, significantly improving the prediction accuracy in sparse data areas.

[0011] Compared to traditional GCN or SA-GCN, which only aggregate local neighborhood information, this model adopts a Graph Transformer architecture and utilizes a multi-head attention mechanism to effectively capture long-range collaborative deformation features within the landslide body. Furthermore, compared to highly complex methods such as ECC or GMMConv that generate huge weight matrices, the SE-PG-GT architecture is more compact and efficient, significantly reducing GPU / memory usage while maintaining high accuracy, and supporting rapid extrapolation of large-scale, high-resolution grid points.

[0012] A highly reliable multi-source data fusion framework was constructed, providing a solid foundation for landslide disaster prevention and control. This invention corrects the planar system bias of InSAR by using high-precision GNSS point references, and combines this with the high-precision reconstruction capabilities of SE-PG-GT to generate a high-precision, physically interpretable planar displacement field. This achievement effectively overcomes the limitations of single monitoring methods, improves the spatial continuity and reliability of landslide monitoring data, and has significant practical implications for accurately identifying landslide hazards and protecting people's lives and property. Attached Figure Description

[0013] Figure 1 This is a schematic diagram of the method flow of the present invention; Figure 2 The annual deformation rate filtering result for the steepest slope in InSAR; Figure 3 A comparison is made between the annual deformation rate residuals of InSAR points optimized by GNSS and the residuals of InSAR points not optimized by GNSS. Figure 4 This is a deformation rate diagram of the Xinpu landslide in 2018. Figure 5 This is the IDW interpolation result for the deformation rate in 2018; Figure 6 The results of SA-GCN interpolation of deformation rates in 2018; Figure 7 The results of the SE-PG-GT model interpolating the deformation rate in 2018 are shown. Figure 8 This is a schematic diagram of the hardware device operation according to an embodiment of the present invention. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0015] Before formally describing the present invention, a general description of the solution of the present invention will be given first to facilitate understanding.

[0016] Example 1 Please refer to Figure 1 The present invention provides a method for interpolating the spatiotemporal characteristics of landslide surface deformation, comprising the following steps: S1. Acquire SAR images of the target landslide body, process them using SBAS-InSAR technology, extract the displacement time series of the landslide body surface, and obtain the original InSAR monitoring point cloud data. S2. The original InSAR monitoring point cloud data is cleaned using the DBSCAN clustering algorithm that integrates spatial and deformation features to identify and remove noise points and outliers, resulting in a cleaned InSAR dataset. It should be noted that step S2 further includes: S21. Take the absolute value of the original InSAR deformation data and construct a significant deformation mask based on its statistical distribution characteristics to screen out potential active areas; Specifically, this invention addresses the original deformation data. Take absolute value Calculate its statistical distribution characteristics and average value. and standard deviation Set a significant deformation threshold. Construct a Boolean mask for significant deformation points .

[0017] As one embodiment, a significant deformation threshold is set. , Alternatively, a fixed value can be set based on historical experience. Construct a Boolean mask for significant deformation points. If a certain point If the number of significant points is 0, it is marked as True; otherwise, it is marked as False. If the number of significant points selected is 0, an empty mask is returned, indicating that the region has no obvious deformation or that the threshold needs to be adjusted.

[0018] S22. Normalize the spatial coordinates and deformation values ​​respectively, and introduce weight factors to construct a feature matrix that integrates spatial proximity and deformation similarity. Specifically, this invention relates to spatial coordinates and deformation value Normalization was performed separately to obtain and Constructing feature vectors ,in This is a weighting factor (e.g., 0.5), used to adjust the contribution ratio of spatial features and deformation features in clustering.

[0019] S23. Apply the DBSCAN algorithm to cluster the fused feature matrix, identify and retain the main deformation clusters that represent the main active area of ​​the landslide, and remove isolated noise points in the clustering results. Specifically, this invention is based on DBSCAN (Density-Based Spatial Clustering of Applications with Noise) for high-deformation cluster identification. The DBSCAN algorithm analyzes the feature matrix... Clustering is performed. Core clusters are identified in the clustering results, and labeled noise points are removed. The cluster with the most points or the highest average deformation rate is defined as the principal deformation cluster, which is the main active area of ​​the landslide. The rest are non-principal deformation areas or background areas.

[0020] S24. For each InSAR point, local anomaly detection is performed within its Euclidean distance neighborhood. Local anomalies are marked by comprehensively evaluating the spatial distance score and deformation anomaly score. Specifically, this invention searches for a set of neighboring points within a certain Euclidean distance range for each InSAR point. Spatial distance score and deformation anomaly score are calculated. If the combined anomaly score of a point exceeds a preset threshold, it is marked as an anomaly and recorded in an anomaly mask. middle.

[0021] S25. Based on the main deformation cluster classification and local anomaly point labeling results, the point cloud is divided into different categories, and median filtering, spatial constraint inverse distance weighting, or moving average smoothing are used for targeted processing to output a cleaned and smoothed InSAR dataset.

[0022] Specifically, based on the main deformation cluster classification in S23 and the anomaly marking in S24, the present invention processes the following three cases respectively: background noise or non-main region anomalies, anomalies in the main deformation region, and normal points in the main deformation region.

[0023] Case 1, background noise or outliers in non-main areas: directly replace with neighborhood median filter to eliminate random noise; Scenario 2, outliers within the main deformation zone: To preserve the deformation characteristics of the main landslide body while eliminating outliers, the spatially constrained inverse distance weighted method is used to reconstruct the value of this point. First, the spatial weights are calculated separately. (Based on the reciprocal of metric distance) and deformation similarity weights (Based on the distribution probability of deformation values ​​at neighboring points). Next, weight fusion is performed: (The weighting ratio can be adjusted based on prior knowledge). Finally, the fusion weights are used to perform a weighted average of the neighboring points, replacing the current outlier.

[0024] Case 3, Normal points within the main deformation region: Slight smoothing is performed using local statistics (such as small window moving average) to preserve the original deformation trend. Output the cleaned InSAR dataset. Obtain the InSAR dataset. Then, calculate the displacement along the steepest slope. The solution equation is as follows: In the formula: It is the velocity of the satellite measured along its line of sight. The slope of the ground. The direction of the surface slope; The incident angle for satellite imaging; The angle between the satellite's orbital direction and geographic north; The angle between the geographic due west direction and the line of sight for satellite imaging; Please refer to Figure 2 , Figure 2 This is the result of filtering the annual deformation rate at the steepest slope in InSAR; S3. Obtain observation data from GNSS stations in the study area, and use KDTree to construct local neighborhood indexes for each GNSS station. Based on the Gaussian kernel function, perform adaptive weighted correction on InSAR points falling into the neighborhood to eliminate systematic bias in InSAR data, realize spatial fusion of multi-source data, and obtain the corrected InSAR point dataset. It should be noted that this step aims to leverage the high-precision location advantage of GNSS data to correct potential biases in InSAR data over long wavelengths, which cover a large area. Through local constraints and smooth transition strategies, the advantages of both data types are complemented.

[0025] It should be noted that step S3 further includes: S31. The time series of the acquired GNSS station observation data is preprocessed to remove gross errors, and downsampled according to the acquisition time point of the InSAR image to obtain GNSS displacement data aligned with the InSAR time series. Specifically, because the sampling intervals of InSAR and GNSS monitoring values ​​differ—InSAR typically uses a 12-day interval, while GNSS uses an hourly interval—effective fusion and comparison of the two data sets requires downsampling the high-frequency GNSS displacement time series to a frequency matching the InSAR data sampling interval. Therefore, GNSS data undergoes downsampling processing. (1) Selecting the nearest InSAR sampling time point: Obtaining the time points of all InSAR observations For each GNSS observation time point Find the closest one in the set of InSAR sampling time points. .

[0026] (2) Extract the GNSS displacement at the corresponding time point: at the found nearest InSAR sampling time point Extract the GNSS station displacement value corresponding to that time point. GNSS displacement time series aligned with InSAR sampling time points were obtained. ,in This represents the number of corresponding points for valid GNSS observations.

[0027] If GNSS data still fluctuates significantly between adjacent InSAR sampling points, a short-term smoothing filter (e.g., a three- or seven-day moving average) can be applied to the downsampled GNSS displacement sequence to reduce short-period noise and make it closer to the slow deformation trend reflected by InSAR. However, over-smoothing may result in the loss of true instantaneous GNSS deformation information.

[0028] S32. Construct the InSAR point data to be corrected into a KDTree spatial index, and quickly retrieve all InSAR points in its local neighborhood by setting a search radius with each GNSS station as the center. As one implementation, the InSAR point data to be corrected is constructed into a KDTree (K-Dimensional Tree) spatial index structure. This is based on the spatial coordinates of each GNSS station. Set the search radius around the center. (This radius is set according to the deformation characteristics of the study area, for example, 10m-200m). KDTree is used to quickly retrieve the search radius falling on each GNSS station. All InSAR points within the range constitute the local InSAR neighborhood set of that GNSS point. .

[0029] S33. Calculate the mean value of InSAR deformation rate in the local neighborhood and compare it with the annual average deformation rate of the central GNSS station. Based on the absolute value of the difference and the preset consistency threshold, determine whether correction is needed and determine the basic deviation amount. As one example, since the single-point measurement accuracy of GNSS is higher than that of InSAR points, the confidence level of the GNSS point monitoring values ​​is set to the highest, and the local InSAR neighborhood set of each GNSS point is calculated. Mean value of InSAR deformation rate calculation Calculate the absolute value of the difference between the annual average deformation rate of GNSS and the local InSAR mean: Set a consistency threshold. (For example, based on the nominal mean square error settings of InSAR and GNSS). Judgment logic: (1) If If it is determined that there is a significant systematic conflict in the InSAR data of the GNSS point or its local area (such as InSAR unwrapping error), the InSAR point in the area is marked, and the correction is carried out in step S34 according to the GNSS monitoring value. (2) If The GNSS point was determined to have a high confidence level, and the baseline deviation was calculated. Proceed to step S34 for correction.

[0030] S34. For the neighborhood of a GNSS point that needs correction, calculate the distance attenuation weight of each InSAR point in its neighborhood based on the Gaussian kernel function, and calculate the single-point correction amount in combination with the basic deviation amount; for InSAR points that fall into the overlapping neighborhood of multiple GNSS points, use the cumulative normalization strategy to perform comprehensive correction and obtain the corrected InSAR deformation rate.

[0031] As one example, calculating InSAR points Euclidean distance to the current GNSS point .

[0032] (1) Construct the Gaussian decay weight function: in, The bandwidth parameter for the Gaussian kernel is set to one-third of the search radius (i.e., This ensures that the correction weight decays smoothly from the center of the GNSS point outwards, avoiding a "step effect" at the correction boundary.

[0033] (2) Calculate the contribution of single-point correction: .

[0034] (3) Multi-GNSS point integrated correction: Considering that an InSAR point may be located in the overlapping area of ​​the search radii of multiple GNSS points, an accumulation normalization strategy is adopted: Accumulation weight: Cumulative correction amount: Calculate the final corrected InSAR deformation rate : (in To prevent tiny quantities with a denominator of zero.

[0035] (4) Output the corrected InSAR point dataset that incorporates GNSS reference information.

[0036] Please refer to Figure 3 , Figure 3 A comparison is made between the annual deformation rate residuals of InSAR points optimized by GNSS and the residuals of InSAR points not optimized by GNSS; from Figure 3 As can be seen, after introducing the GNSS reference constraint, the mean residual decreased from 4.38 mm to 4.13 mm, and the main peak of the error distribution converged significantly towards zero. This characteristic indicates that this step effectively reduced the systematic bias in InSAR observations and improved the overall absolute accuracy of deformation monitoring.

[0037] S4. Construct a heterogeneous fully connected graph containing geometric coordinates and multi-source environmental factors, and use the corrected InSAR point data and the grid points to be interpolated as graph nodes to form graph structure data. It should be noted that this step aims to utilize deep graph neural network technology to establish a nonlinear mapping relationship from sparse, discrete InSAR observation points to a global continuous raster deformation field. This model integrates the long-range spatial dependency capture capability of the Graph Transformer, the feature adaptive filtering capability of the compression-excitation (SE) module, and prior knowledge of terrain physical constraints.

[0038] Step S4 further includes: S41. For all nodes within the study area, including known InSAR points and grid points to be interpolated, extract elevation, slope, and aspect environmental factors from the digital elevation model (DEM), and encode periodic aspect data using trigonometric functions, then concatenate them with normalized geometric coordinates to construct high-dimensional node feature vectors. As one example, multi-source environmental feature engineering and node multi-dimensional vector construction are necessary to build a high-dimensional feature vector that includes geometric location and geological environmental attributes in order to enable the model to recognize the landslide-prone environment.

[0039] (1) Environmental factor extraction: For each node in the study area (include Training nodes with known deformation and (one prediction node to be interpolated), obtain its two-dimensional plane coordinates. And extract the corresponding elevation from the digital elevation model (DEM). ,slope and slope (Unit: degrees, range) ).

[0040] (2) Slope aspect periodicity coding: Given the circular periodicity of slope aspect data (e.g. and The only difference in physical space However, the numerical differences are huge. To ensure numerical continuity, trigonometric function encoding is used to decompose the slope aspect into two orthogonal components: (3) Feature normalization and splicing: for coordinates Perform Min-Max normalization and map to Interval. For environmental factors. Standardize it so that its mean is 0 and its variance is 1.

[0041] By concatenating the above features, a dimension of [value missing] is formed. The node input feature vector : .

[0042] S42. Using all nodes as vertices, based on the Euclidean distance between nodes, use the K-nearest neighbor algorithm to search for and connect the K nearest neighbors of each node, and construct a heterogeneous KNN graph structure that includes node features and edge geometric attributes.

[0043] Specifically, a heterogeneous KNN graph containing geometric properties is constructed to simulate the material transport and force transmission within the landslide body, and a directed graph structure is also built. First, a fully connected graph is constructed. Then, the K-nearest neighbor algorithm is used, with the Euclidean distance between nodes as the metric, to search for the nearest neighbor for each node. Neighbors (in the example) Edge attribute extraction: For each edge from the source node... Point to target node edge Calculate its geometric attribute vector To capture relative positional relationships: in Let be the Euclidean distance between the two points. This vector is then normalized to its maximum value to prevent excessively large values ​​from causing gradient explosion.

[0044] S5. Construct and train the Physical Guided Graph Transformation Network (SE-PG-GT) model based on channel attention enhancement. Input the graph structure data into the SE-PG-GT model and learn the mapping relationship from sparse discrete points to continuous grids under the semi-supervised transductive learning framework. The SE-PG-GT model in step S5 consists of multiple stacked GT-SE modules, and the function of each module includes: First, a multi-head self-attention mechanism is used to aggregate information from neighboring nodes, where the calculation of the attention coefficient considers both node features and edge geometric properties to capture long-distance spatial dependencies and nonlinear relationships between nodes. Then, the aggregated features are input into the compression-excitation SE module. Through global feature compression and channel weight excitation and recalibration, the feature channels that play a key role in landslide deformation are dynamically selected and enhanced.

[0045] Specifically, the SE-PG-GT network architecture is constructed, with the main network component consisting of... It consists of stacked GT-SE Blocks (Graph Transform-Channel Attention Modules). In this embodiment... The specific hierarchical operations are as follows: (1) Spatial Attention Feature Aggregation: The first layer utilizes a multi-head attention mechanism to capture the spatial dependencies between nodes. Let the input feature dimension be... The output dimension is The number of attention heads is (In the example) For nodes and his neighbors First calculate the first A query for each attention head ( ), Key( ) and Value( Vector: in This is a learnable weight matrix. Attention coefficients are calculated using a dot-product scaling model. : node aggregation features Piecing together weighted information for all neighbors: at this time, The dimension is Then BatchNormalization and ELU activation are performed.

[0046] (2) Channel Attention Calibration (SE-Block): To enable the model to automatically identify key environmental factors (such as slope features being more important in steep areas), a compression-excitation mechanism is introduced to reweight the feature channels. The following is a detailed implementation of compression, excitation, and reweighting in this example: Compression, utilizing fully connected layers From feature dimensions Compress to ( ), and activated by ReLU: Incentives, utilizing fully connected layers The dimensions are restored, and channel weight vectors are generated using the Sigmoid function. (scope ): Reweighting will generate a weight vector. With original features Perform element-wise product: (3) Network layer stacking Layer 1: Input 6 dimensions, output Dimension, extract low-level geometric features.

[0047] Layer 2: Input 256 dimensions, output Dimensions enable deep semantic interaction.

[0048] Layer 3: Input 256 dimensions, output Dimension (single-head attention), fusion features.

[0049] Output layer: Linear mapping layer, which maps 64-dimensional features to 1-dimensional predicted deformation rate. .

[0050] It should be noted that the loss function used in step S5 to train the SE-PG-GT model is a hybrid loss function that includes physical consistency constraints, and its construction method is as follows: Construct a data fitting loss term and calculate the mean square error between the model's predicted values ​​and the actual values ​​based only on known GNSS / InSAR points; Specifically, data fitting loss ( ): Only for the training node set (Given points) Calculate the predicted value Compared with the true value Mean square error: .

[0051] A terrain physical constraint loss term is constructed. Based on the principle of landslide dynamics, a penalty function is constructed using the normalized slope value. The predicted deformation values ​​of all nodes in the entire map, including unknown grid points, are constrained, and large deformation predictions generated in flat areas are penalized. Specifically, terrain physical constraint loss ( Based on the principle of landslide dynamics—"flat terrain makes it difficult to accumulate gravitational potential energy, thus preventing severe deformation"—this constraint applies to the entire map node set. (Including unknown points). Construct the penalty function: in The normalized slope value ( ), This is the constraint strength coefficient (3.0 in the example).

[0052] When the slope (Steep) The penalty is minimal, allowing the model to predict large deformations.

[0053] When the slope When (flat), The penalty is extremely high, forcing the model's predicted values ​​to be... Approaching 0.

[0054] The model is trained by adding the data fitting loss term and the terrain physical constraint loss term using a balancing coefficient to form the final overall optimization objective.

[0055] Specifically, the overall optimization objective is: in This is the balance coefficient.

[0056] S6. Use the trained SE-PG-GT model to predict all nodes in the map and output a high-precision, physically constrained, continuous planar deformation field covering the entire landslide body.

[0057] It should be noted that step S6 further includes: With the parameters of the trained SE-PG-GT model fixed, forward propagation calculation is performed on all nodes in the full graph containing all grid points to be interpolated. Extract the grid point prediction values ​​from the model output and perform inverse standardization to obtain the deformation rate in real physical units; The one-dimensional prediction results are reshaped into a two-dimensional matrix according to the preset grid coordinates to generate a high-resolution displacement rate distribution map covering the entire landslide body.

[0058] Specifically, the entire map data is input into the network, and the Adam optimizer is used for iterative parameter updates until the total loss function converges. Keeping the network parameters fixed, the network's output values ​​for the prediction nodes (grid points) are extracted. Finally, the output values ​​are output and inversely standardized (Inverse StandardScaler) to obtain the true deformation rate in mm / year. The one-dimensional prediction results are then reconstructed into a two-dimensional matrix according to the grid coordinates, generating a high-resolution displacement rate distribution map covering the entire landslide body.

[0059] Following step S6, a model evaluation step is also included: Select preset test points, and calculate the root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²). 2 The accuracy of the generated high-resolution displacement rate distribution map is quantitatively evaluated.

[0060] Specifically, the root mean square error (RMSE) is calculated using the following formula: ; Calculate the mean absolute error (MAE) using the following formula: ); Calculate the coefficient of determination The formula is: .

[0061] Please refer to Figures 4-7 . Figure 4 This is a deformation rate diagram of the Xinpu landslide in 2018. Figure 5 This is the IDW interpolation result for the deformation rate in 2018; Figure 6 The results of SA-GCN interpolation of deformation rates in 2018; Figure 7 The results of the SE-PG-GT model interpolating the deformation rate in 2018 are shown. from Figures 5-7 It can be seen that, compared with the "bull's-eye artifacts" of traditional IDW and the "oversmoothing and boundary leakage" of SAGCN, the proposed SE-PG-GT method overcomes the blindness of two-dimensional space. While accurately capturing extreme deformation values, it restores the true physical boundary of the landslide deformation area constrained by micro-topography. This shows that when dealing with highly nonlinear geological hazards, smoothing algorithms relying solely on spatial distance are ineffective. Real three-dimensional terrain constraints and graph network topology can be introduced to generate reliable deformation surfaces supported by geoscientific mechanisms.

[0062] Example 2: Please see Figure 8 , Figure 8 This is a schematic diagram of the hardware device in operation according to an embodiment of the present invention. The hardware device specifically includes: a landslide surface deformation spatiotemporal feature interpolation device 401, a processor 402, and a storage medium 403.

[0063] A landslide surface deformation spatiotemporal feature interpolation device 401: The landslide surface deformation spatiotemporal feature interpolation device 401 implements the landslide surface deformation spatiotemporal feature interpolation method.

[0064] Processor 402: The processor 402 loads and executes the instructions and data in the storage medium 403 to implement the landslide surface deformation spatiotemporal feature interpolation method.

[0065] Storage medium 403: The storage medium 403 stores instructions and data; the storage medium 403 is used to implement the spatiotemporal feature interpolation method for landslide surface deformation.

[0066] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A landslide surface deformation spatio-temporal feature interpolation method, characterized in that: Includes the following steps: S1. Acquire SAR images of the target landslide body, process them using SBAS-InSAR technology, extract the displacement time series of the landslide body surface, and obtain the original InSAR monitoring point cloud data. S2. The original InSAR monitoring point cloud data is cleaned using the DBSCAN clustering algorithm that integrates spatial and deformation features to identify and remove noise points and outliers, resulting in a cleaned InSAR dataset. S3. Obtain observation data from GNSS stations in the study area, and use KDTree to construct local neighborhood indexes for each GNSS station. Based on the Gaussian kernel function, perform adaptive weighted correction on InSAR points falling into the neighborhood to eliminate systematic bias in InSAR data, realize spatial fusion of multi-source data, and obtain the corrected InSAR point dataset. S4. Construct a heterogeneous fully connected graph containing geometric coordinates and multi-source environmental factors, and use the corrected InSAR point data and the grid points to be interpolated as graph nodes to form graph structure data. S5. Construct and train the Physical Guided Graph Transformation Network (SE-PG-GT) model based on channel attention enhancement. Input the graph structure data into the SE-PG-GT model and learn the mapping relationship from sparse discrete points to continuous grids under the semi-supervised transductive learning framework. S6. Use the trained SE-PG-GT model to predict all nodes in the map and output a high-precision, physically constrained, continuous planar deformation field covering the entire landslide body.

2. The landslide surface deformation spatio-temporal feature interpolation method according to claim 1, characterized in that: Step S2 further includes: S21. Take the absolute value of the original InSAR deformation data and construct a significant deformation mask based on its statistical distribution characteristics to screen out potential active areas; S22. Normalize the spatial coordinates and deformation values ​​respectively, and introduce weight factors to construct a feature matrix that integrates spatial proximity and deformation similarity. S23. Apply the DBSCAN algorithm to cluster the fused feature matrix, identify and retain the main deformation clusters that represent the main active area of ​​the landslide, and remove isolated noise points in the clustering results. S24. For each InSAR point, local anomaly detection is performed within its Euclidean distance neighborhood. Local anomalies are marked by comprehensively evaluating the spatial distance score and deformation anomaly score. S25. Based on the main deformation cluster classification and local anomaly point labeling results, the point cloud is divided into different categories, and median filtering, spatial constraint inverse distance weighting, or moving average smoothing are used for targeted processing to output a cleaned and smoothed InSAR dataset.

3. The landslide surface deformation spatiotemporal feature interpolation method according to claim 1, characterized in that: Step S3 further includes: S31. The time series of the acquired GNSS station observation data is preprocessed to remove gross errors, and downsampled according to the acquisition time point of the InSAR image to obtain GNSS displacement data aligned with the InSAR time series. S32. Construct the InSAR point data to be corrected into a KDTree spatial index, and quickly retrieve all InSAR points in its local neighborhood by setting a search radius with each GNSS station as the center. S33. Calculate the mean value of InSAR deformation rate in the local neighborhood and compare it with the annual average deformation rate of the central GNSS station. Based on the absolute value of the difference and the preset consistency threshold, determine whether correction is needed and determine the basic deviation amount. S34. For the neighborhood of a GNSS point that needs correction, calculate the distance attenuation weight of each InSAR point in its neighborhood based on the Gaussian kernel function, and calculate the single-point correction amount in combination with the basic deviation amount; for InSAR points that fall into the overlapping neighborhood of multiple GNSS points, use the cumulative normalization strategy to perform comprehensive correction and obtain the corrected InSAR deformation rate.

4. The method for interpolating the spatiotemporal characteristics of landslide surface deformation according to claim 1, characterized in that: Step S4 further includes: S41. For all nodes within the study area, including known InSAR points and grid points to be interpolated, extract elevation, slope, and aspect environmental factors from the digital elevation model (DEM), and encode periodic aspect data using trigonometric functions, then concatenate them with normalized geometric coordinates to construct high-dimensional node feature vectors. S42. Using all nodes as vertices, based on the Euclidean distance between nodes, use the K-nearest neighbor algorithm to search for and connect the K nearest neighbors of each node, and construct a heterogeneous KNN graph structure that includes node features and edge geometric attributes.

5. The landslide surface deformation spatiotemporal feature interpolation method according to claim 1, characterized in that: The SE-PG-GT model in step S5 consists of multiple stacked GT-SE modules, and the function of each module includes: First, a multi-head self-attention mechanism is used to aggregate information from neighboring nodes, where the calculation of the attention coefficient considers both node features and edge geometric properties to capture long-distance spatial dependencies and nonlinear relationships between nodes. Then, the aggregated features are input into the compression-excitation SE module. Through global feature compression and channel weight excitation and recalibration, the feature channels that play a key role in landslide deformation are dynamically selected and enhanced.

6. The method for interpolating the spatiotemporal characteristics of landslide surface deformation according to claim 5, characterized in that, In step S5, the loss function used to train the SE-PG-GT model is a hybrid loss function that includes physical consistency constraints, and its construction method is as follows: Construct a data fitting loss term and calculate the mean square error between the model's predicted values ​​and the actual values ​​based only on known GNSS / InSAR points; A terrain physical constraint loss term is constructed. Based on the principle of landslide dynamics, a penalty function is constructed using the normalized slope value. The predicted deformation values ​​of all nodes in the entire map, including unknown grid points, are constrained, and large deformation predictions generated in flat areas are penalized. The model is trained by adding the data fitting loss term and the terrain physical constraint loss term using a balancing coefficient to form the final overall optimization objective.

7. The landslide surface deformation spatiotemporal feature interpolation method according to claim 1, characterized in that: Step S6 further includes: With the parameters of the trained SE-PG-GT model fixed, forward propagation calculation is performed on all nodes in the full graph containing all grid points to be interpolated. Extract the grid point prediction values ​​from the model output and perform inverse standardization to obtain the deformation rate in real physical units; The one-dimensional prediction results are reshaped into a two-dimensional matrix according to the preset grid coordinates to generate a high-resolution displacement rate distribution map covering the entire landslide body.

8. The landslide surface deformation spatiotemporal feature interpolation method according to claim 1, characterized in that: Following step S6, a model evaluation step is also included: Select the preset test point, through the calculation of root mean square error RMSE, mean absolute error MAE and determination coefficient R 2 The accuracy of the generated high-resolution displacement rate distribution map is quantitatively evaluated.

9. A storage medium, characterized in that: The storage medium stores instructions and data to implement the landslide surface deformation spatiotemporal feature interpolation method according to any one of claims 1 to 8.

10. A device for interpolating the spatiotemporal characteristics of landslide surface deformation, characterized in that: include: A processor and a storage medium; the processor loads and executes instructions and data in the storage medium to implement the landslide surface deformation spatiotemporal feature interpolation method according to any one of claims 1 to 8.

Images