Intelligent remote sensing interpretation combined with topographic factor geological disaster hidden danger point accurate positioning device

Abstract

This invention discloses a device for precise location of geological hazard hazard points by combining intelligent remote sensing interpretation with topographic factors. It relates to the field of geological hazard location technology and includes: a data acquisition terminal for acquiring multi-temporal synthetic aperture radar interferometric image sequences and digital elevation model data, performing spatiotemporal baseline registration to generate a temporal interferometric phase dataset, and calculating a set of topographic factors; a data processing center for separating weak deformation signals through variational mode decomposition and independent component analysis, fusing denoised deformation phase sequences and topographic factors to generate a multi-dimensional spatiotemporal feature tensor, constructing a topographic topology map, inputting it into a graph attention network model to output hazard probability values ​​and boundary confidence, determining the location of hazard points, and tracking the boundaries of hazard bodies; and a remote access terminal for performing visualization and spatial database storage. This invention effectively improves the early and accurate identification capability of large, deep landslide hazards with strong concealment and the accurate delineation capability of hazard body boundaries.

Description

Technical Field

[0001] This invention relates to the field of geological disaster location technology, and in particular to a device for precise location of geological disaster hazard points by combining intelligent remote sensing interpretation with topographic factors. Background Technology

[0002] Precise location of potential geological hazard sites is a crucial prerequisite for disaster prevention and mitigation. Currently, the accurate location of potential large-scale deep landslide hazards mainly relies on a collaborative technology system integrating air, space, and ground. This system typically follows a tiered technical approach: first, wide-area scanning is conducted using spaceborne synthetic aperture radar interferometry to identify areas of abnormal surface deformation; then, technologies such as UAV-borne lidar are used to conduct detailed investigations of the anomalous areas, obtaining high-precision topographic data; finally, verification and continuous monitoring are achieved through the deployment of ground-based and borehole monitoring equipment.

[0003] However, when dealing with potential large-scale deep landslides, these landslides often exhibit large-scale, uniform, slow-moving, and homogeneous deep creep in their early stages, with extremely weak surface deformation signals. Existing wide-area scanning methods based on spaceborne synthetic aperture radar often detect surface deformation levels at the same level as system noise and atmospheric delay errors. This makes it extremely difficult to effectively separate and confirm this spatially continuous, weak, true deformation pattern from complex interferometric signals, easily leading to missed or false diagnoses. Furthermore, the macroscopic topographic features of such landslides are usually insignificant, and may only appear as gentle slopes on conventional high-precision digital elevation models. This reduces the effectiveness of subsequent topographic feature-based identification models due to a lack of effective input features, making it difficult to accurately delineate the complete boundaries of the hazard body.

[0004] Therefore, existing technical solutions face core technical challenges in addressing potential large-scale deep landslide hazards. These challenges include the difficulty in reliably extracting weak early deformation signals from noise, and the weak indicative power of surface topography, leading to insufficient early identification accuracy and inaccurate hazard boundary delineation. A new technical solution is urgently needed to improve the early and accurate identification capability of such highly concealed and hazardous geological hazards. Summary of the Invention

[0005] The purpose of this invention is to provide a device for accurately locating geological hazard points by combining intelligent remote sensing interpretation with topographic factors, thereby solving the problems of the prior art.

[0006] To achieve the above objectives, the present invention provides a device for precise location of geological hazard points by combining intelligent remote sensing interpretation with topographic factors, comprising: a data acquisition terminal, a data processing center, and a remote access terminal; The data acquisition terminal is used to perform the following steps: Acquire multi-temporal synthetic aperture radar (SAR) interferometric image sequences of the target area, perform spatiotemporal baseline registration on the multi-temporal SAR SAR interferometric image sequences, and generate a temporal interferometric phase dataset; acquire digital elevation model (DEM) data of the target area, and calculate a set of terrain factors based on the DEM data; The data processing center is used to perform the following steps: The system receives a temporal interferometric phase dataset and a set of topographic factors. It performs variational mode decomposition (VMD) on the temporal interferometric phase dataset pixel by pixel, decomposing the temporal interferometric phase of each pixel into multiple intrinsic mode components (IMCs). An observation matrix is ​​constructed from these IMCs, and independent component analysis (ICA) is performed on the observation matrix to separate the atmospheric delay IMC, system noise IMC, and deformation signal IMC. The deformation signal IMC is extracted and reconstructed to generate a denoised deformation phase sequence. Spatial registration is performed between the denoised deformation phase sequence and the set of topographic factors to generate a multidimensional spatiotemporal feature tensor. A topographic topology map is constructed based on the spatial grid division of the target area, and the multidimensional spatiotemporal feature tensor is mapped to the nodes of the topographic topology map to generate an attributed topographic feature map. The attributed topographic feature map is input into a pre-trained graph attention network (GAN) model, which outputs the hazard probability value and boundary confidence of each node. Nodes with hazard probability values ​​greater than a preset probability threshold are selected to determine the spatial location of geological hazard points. Regional growth boundary tracking is performed on adjacent hazard nodes based on the boundary confidence to generate hazard body boundary vector data. The remote access terminal is used to perform the following steps: Receive spatial location data of geological hazard hazard points and vector data of hazard body boundaries; overlay the spatial location data of geological hazard hazard points and vector data of hazard body boundaries onto the geographic information system base map to generate a hazard distribution visualization layer; store the hazard distribution visualization layer into a spatial database.

[0007] Furthermore, the generation of the temporal interferometric phase dataset includes: selecting the main image from the multi-temporal synthetic aperture radar interferometric image sequence and using the remaining images as auxiliary images; calculating the offset between images using the intensity cross-correlation method, establishing the geometric mapping relationship between the auxiliary images and the main image through polynomial fitting, resampling the auxiliary images into the geometric framework of the main image to eliminate the geometric offset between images; calculating the temporal baseline and spatial baseline between each auxiliary image and the main image, where the temporal baseline is defined as the difference between the acquisition time of the auxiliary image and the acquisition time of the main image, and the spatial baseline is defined as the distance between the satellite orbit position of the auxiliary image and the main image. The satellite orbital position is projected at a distance perpendicular to the radar line of sight; image pairs that meet the criteria of a spatial baseline less than two-thirds of a preset critical baseline and a temporal baseline less than a preset time decoherence threshold are selected to construct an interferometric pair set; the complex values ​​of the main image and the auxiliary image are multiplied by their conjugates to obtain an interferometric complex map; the phase information of the interferometric complex map is extracted to obtain the entangled interferometric phase; differential interferometry is used to sequentially remove the flat land phase and the terrain phase from the original interferometric phase to obtain a differential interferometric phase map; the interferometric phase maps are arranged in chronological order to generate a time-series interferometric phase dataset.

[0008] Furthermore, the set of terrain factors is calculated based on the digital elevation model data, including: performing a sliding window traversal on the digital elevation model data to calculate the slope and aspect data for each pixel, where the slope data is the rate of change of elevation and the aspect data is the azimuth of the direction of maximum slope; calculating the planar curvature and profile curvature data for each pixel based on the second-order partial derivatives of the digital elevation model data; calculating the difference between the maximum and minimum elevations within a preset neighborhood window to obtain terrain relief data; calculating the standard deviation of the elevation values ​​within the neighborhood window to obtain terrain roughness data; and combining the slope data, aspect data, planar curvature data, profile curvature data, terrain relief data, and terrain roughness data into a set of terrain factors.

[0009] Furthermore, performing variational mode decomposition on a pixel-by-pixel basis on the temporal interferometric phase dataset includes: constructing a constrained optimization problem for variational mode decomposition, taking the phase time series of each pixel in the temporal interferometric phase dataset as the signal to be decomposed; iteratively solving the constrained optimization problem using the alternating direction multiplier method to obtain the center frequency and bandwidth of each intrinsic mode component; and performing adaptive decomposition on the phase time series based on the center frequency and bandwidth to output K intrinsic mode components, where K is a positive integer.

[0010] Furthermore, performing independent component analysis on the observation matrix includes: constructing an observation matrix by arranging the K intrinsic mode components in rows; performing centering and whitening processes on the observation matrix sequentially to obtain a whitened matrix; performing independent component analysis on the whitened matrix based on the non-Gaussianity maximization criterion to obtain an unmixed matrix; separating the observation matrix into statistically independent components through the unmixed matrix; and identifying each independent component as an atmospheric delay independent component, a system noise independent component, and a deformation signal independent component based on their spatiotemporal characteristics.

[0011] Furthermore, the generation of the multidimensional spatiotemporal feature tensor includes: uniformly resampling the denoised deformation phase sequence and the set of terrain factors to the same spatial resolution; expanding the resampled denoised deformation phase sequence along the time dimension to obtain the deformation feature matrix; stacking the terrain factor data of each terrain factor in the set along the feature channel dimension to obtain the terrain feature matrix; and concatenating the deformation feature matrix and the terrain feature matrix along the feature channel dimension to generate the multidimensional spatiotemporal feature tensor.

[0012] Furthermore, the generation of the attributed terrain feature map includes: dividing the target area into regular spatial grid cells; defining a set of nodes for the terrain topology map, with each node corresponding to a spatial grid cell; defining a set of edges for the terrain topology map, with adjacent grid cells connected by an eight-neighborhood connection rule; and assigning the feature vector corresponding to each grid cell in the multidimensional spatiotemporal feature tensor to the corresponding node of the terrain topology map to generate the attributed terrain feature map.

[0013] Furthermore, the attributed terrain feature map is input into a pre-trained graph attention network model, which outputs the hazard probability value and boundary confidence of each node. The graph attention network model includes an input layer, multiple graph attention convolutional layers, and an output layer. The input layer receives the node feature vectors from the attributed terrain feature map. Each graph attention convolutional layer calculates the attention weight coefficients between the current node and its neighboring nodes through a multi-head attention mechanism, performs weighted aggregation on the features of neighboring nodes based on the attention weight coefficients, and updates the hidden state representation of the current node. The output layer performs classification mapping on the hidden state representation of the nodes output by the last graph attention convolutional layer through two parallel fully connected branches, and outputs the hazard probability value and boundary confidence of each node.

[0014] Furthermore, based on the boundary confidence level, the region growth boundary tracking is performed on adjacent hazard nodes to generate hazard body boundary vector data, including: selecting the hazard candidate node with the highest hazard probability value as the seed node; starting from the seed node, judging whether the neighboring nodes belong to the same hazard body based on the boundary confidence level, and including the neighboring nodes that meet the boundary confidence level condition into the current hazard body region; iteratively performing region growth until the boundary of the current hazard body region converges; extracting the boundary node sequence of the outer contour of the current hazard body region, converting the coordinate sequence of the boundary nodes into vector polygon format, and generating hazard body boundary vector data.

[0015] Furthermore, the criteria for determining whether a neighboring node belongs to the same hazard body based on the boundary confidence level are as follows: if the boundary confidence level of a neighboring node is greater than the preset boundary confidence level threshold, and the neighboring node belongs to the set of hazard candidate nodes, then the neighboring node belongs to the same hazard body; otherwise, the neighboring node does not belong to the same hazard body.

[0016] Therefore, the present invention employs the above-mentioned intelligent remote sensing interpretation combined with terrain factors for precise location of geological disaster hazard points, which has the following beneficial effects: 1. By performing variational mode decomposition on a pixel-by-pixel basis on the temporal interferometric phase dataset, the complex interferometric phase signal is decomposed into multiple intrinsic mode components with different center frequencies. Then, independent component analysis is performed on multiple intrinsic mode components. Based on the statistical independence criterion, atmospheric delay error, system noise and real deformation signal are separated into independent components. This enables reliable extraction of weak deformation signals in the early stage of large deep landslides and solves the technical problem that weak deformation signals and noise are of similar magnitude and difficult to separate effectively in the existing technology. 2. By fusing the denoised deformation phase sequence with the set of terrain factors to generate a multidimensional spatiotemporal feature tensor, and constructing a terrain topology map input map attention network model, the spatial dependency relationship between each grid unit is adaptively learned using a multi-head attention mechanism. This can uncover hidden terrain patterns with spatial continuity but indistinct local features, solving the technical problem of weak surface terrain feature indicativeness leading to reduced recognition model performance in existing technologies. 3. By combining variational mode decomposition and independent component analysis, high-quality denoised input features are provided for the graph attention network model. The spatial correlation learning ability of the graph attention network model further explores the deep correlation between the denoised deformation signal and multidimensional terrain factors, which improves the early and accurate identification ability of large-scale deep landslide hazards with strong concealment and the accurate delineation of the boundary of the hazard body.

[0017] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0018] Fig. 1This is a flowchart illustrating the execution process of the intelligent remote sensing interpretation combined with terrain factors for precise location of geological hazard points according to the present invention. Fig. 2 This is a schematic diagram of the structure of the intelligent remote sensing interpretation combined with topographic factors for precise location of geological disaster hazard points according to the present invention. Detailed Implementation

[0019] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely illustrates selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0020] Please see Figs. 1-2 The device for accurately locating geological hazard points by combining intelligent remote sensing interpretation with topographic factors includes: a data acquisition terminal 1, a data processing center 2, and a remote access terminal 3.

[0021] In an embodiment of the present invention, the monitoring needs of a potential large-scale deep landslide hazard area in a mountainous region are addressed. This area covers approximately 50 square kilometers, with significant topographic relief and slopes ranging from 5 to 45 degrees, and multiple suspected signs of deep creep exist. Traditional spaceborne synthetic aperture radar interferometry methods struggle to effectively identify early, weak deformation signals in this area. The technical solution provided by the present invention effectively addresses this problem.

[0022] Data acquisition terminal 1 is configured to acquire multi-temporal synthetic aperture radar interferometric image sequences of the target area, perform spatiotemporal baseline registration on the multi-temporal synthetic aperture radar interferometric image sequences to generate a temporal interferometric phase dataset; acquire digital elevation model data of the target area, and calculate a set of terrain factors based on the digital elevation model data. The specific execution process includes the following steps: Step 101: Acquisition of multi-temporal synthetic aperture radar interferometric image sequence and spatiotemporal baseline registration; Data acquisition terminal 1 acquires a multi-temporal synthetic aperture radar (SAR) interferometric image sequence of the target area from the satellite data distribution center. In an embodiment of the present invention, the acquired image sequence includes 36 C-band SAR single-look complex images, spanning 36 months, with a revisit period of 12 days.

[0023] The specific process of performing spatiotemporal baseline registration on multi-temporal synthetic aperture radar interferometric image sequences is as follows: Step A1: Select the main image from the multi-temporal synthetic aperture radar interferometric image sequence. The selection criteria are: select the image located in the middle of the sequence with good atmospheric conditions as the main image, and the remaining images as auxiliary images.

[0024] Step A2: Perform spatial registration between each auxiliary image and the main image. The offset between images is calculated using the intensity cross-correlation method. A geometric mapping relationship between the auxiliary and main images is established through polynomial fitting. The auxiliary images are resampled to fit within the geometric framework of the main image, eliminating the geometric offset between images. The registration accuracy is controlled within 0.1 pixels.

[0025] Step A3: Calculate the temporal and spatial baselines between each secondary image and the primary image. Temporal baseline Defined as the difference between the acquisition time of the secondary image and the acquisition time of the primary image, in days. Spatial baseline Defined as the projected distance between the orbital positions of the secondary image satellite and the primary image satellite in a direction perpendicular to the radar line of sight, in meters.

[0026] Step A4: Select image pairs whose temporal and spatial baselines satisfy the coherence condition, and construct an interferometric pair set. The coherence condition is: spatial baseline... Less than two-thirds of the critical baseline, time baseline The time decoherence threshold is less than the threshold value. In an embodiment of the present invention, the spatial baseline threshold is set to 200 meters and the time baseline threshold is set to 180 days, resulting in 215 valid interference pairs after screening.

[0027] Step A5: Perform interferometric processing on each interferometric pair in the interferometric pair set to generate an interferometric phase map. The interferometric processing includes: multiplying the complex values ​​of the main image and the auxiliary image by their conjugates to obtain a complex interferometric map; extracting the phase information from the complex interferometric map to obtain the wrapped interferometric phase; and removing the flat terrain phase and topographic phase to obtain a differential interferometric phase map. Differential Interferometric Phase The calculation formula is: ; in, The original interference phase, It is a flat phase. This refers to the terrain phase.

[0028] Step A6: Arrange the interferometric phase images in chronological order to generate a time-series interferometric phase dataset. The dataset has a three-dimensional array structure with dimensions of time (number of interferometric pairs), row (number of image rows), and column (number of image columns).

[0029] Step 102, Calculation of terrain factor set; Data acquisition terminal 1 acquires digital elevation model data of the target area, which may be derived from measurements by aerospace mapping satellites or airborne lidar. In this embodiment of the invention, digital elevation model data with a resolution of 12.5 meters is used.

[0030] The specific process for calculating the set of terrain factors based on digital elevation model data is as follows: Step B1: Calculate slope data. Perform a 3×3 sliding window traversal on the digital elevation model data, calculating the slope value of the center pixel within each window. Assume the elevation values ​​of the 9 pixels within the window are arranged in order of position as follows: to ,in The central pixel has a pixel spacing of . (Unit: meters). Slope The calculation formula is: ; ; ; in, It is the arctangent function. The elevation change rate in the east-west direction. This represents the rate of change of elevation in the north-south direction.

[0031] Step B2: Calculate slope aspect data. Slope aspect The azimuth angle for the direction of maximum gradient is calculated by taking true north as 0 degrees and increasing clockwise. The formula is as follows: ; when hour, ;when and hour, .

[0032] Step B3: Calculate planar curvature data and profile curvature data. Calculate curvature based on the second-order partial derivatives of the digital elevation model data. Let... , , The elevation function is respectively in Second-order partial derivatives of the direction, The second-order partial derivatives and mixed partial derivatives in the direction are calculated using the following formulas: ; ; ; Plane curvature The formula for representing the curvature of contour lines is: ; Section curvature The formula for representing the curvature change along the slope direction is: ; Step B4: Calculate terrain relief data. Calculate the maximum elevation within a preset neighborhood window (a 21×21 pixel window is used in this embodiment). With minimum elevation The difference is used to obtain the terrain relief. .

[0033] Step B5: Calculate terrain roughness data. Calculate the standard deviation of elevation values ​​within the neighborhood window. As terrain roughness : ; in, This represents the total number of pixels within the window. This represents the average elevation within the window.

[0034] Step B6: Combine the slope data, aspect data, plane curvature data, profile curvature data, terrain relief data, and terrain roughness data into a terrain factor set. The data structure of this set is a three-dimensional array, with dimensions of feature channel dimension (6 terrain factors), row dimension, and column dimension.

[0035] Data processing center 2 is configured to receive the time-series interferometric phase dataset and terrain factor set transmitted by data acquisition terminal 1, and perform signal decomposition, feature fusion and intelligent recognition processing. The specific execution process includes the following steps: Step 201, variational mode decomposition; Variational mode decomposition (VMD) is performed pixel-by-pixel on the temporal interferometric phase dataset, decomposing the temporal interferometric phase of each pixel into multiple intrinsic mode components. VMD is an adaptive signal decomposition method that can decompose complex signals into narrowband components with different center frequencies.

[0036] Assuming the temporal interferometric phase sequence of each pixel is taken as the original signal, it needs to be decomposed into One intrinsic mode component It is a positive integer, and in this embodiment, it is taken as... Variational mode decomposition (VMD) constructs a constrained optimization problem to minimize the sum of estimated bandwidths of all intrinsic mode components (IMCs) while ensuring that the sum of all IMCs equals the original signal. By introducing Lagrange multipliers and a quadratic penalty factor, the constrained optimization problem is transformed into an augmented Lagrange function. An alternating direction multiplier method is used to iteratively update the IMCs, center frequency, and Lagrange multipliers in the frequency domain until the difference between two adjacent iterations is less than a preset convergence threshold.

[0037] In an embodiment of the present invention, variational mode decomposition is performed on the 215 temporal interference phase values ​​of each pixel, with a penalty factor... Set to 2000, update step size Set to 0, and the convergence threshold to [value missing]. The decomposition yields five intrinsic mode components, each corresponding to a signal component with different frequency characteristics.

[0038] Step 202, Independent component analysis; An observation matrix is ​​constructed for multiple intrinsic mode components. Independent component analysis is then performed on the observation matrix to separate the atmospheric delay independent component, the system noise independent component, and the deformation signal independent component. The specific steps include: Step C1: The intrinsic mode components are arranged in rows to construct the observation matrix, with dimension . ,in The number of time sampling points included for each intrinsic mode component.

[0039] Step C2: Perform centering on the observation matrix, specifically by calculating the mean of each row and subtracting it to obtain a zero-mean observation matrix.

[0040] Step C3: Perform whitening on the zero-mean observation matrix. Specifically, the eigenvalue decomposition method is used to calculate the eigenvectors and eigenvalues ​​of the covariance matrix of the zero-mean observation matrix, constructing a whitening transformation matrix. The zero-mean observation matrix is ​​then whitened using this whitening transformation matrix, making the components uncorrelated and normalizing the variance, thus obtaining the whitened matrix.

[0041] Step C4: Using the FastICA algorithm, independent component analysis is performed on the whitening matrix based on the non-Gaussianity maximization criterion. Negative entropy is used as the non-Gaussianity metric, and a fixed-point iterative algorithm is employed to solve for the unmixing vectors corresponding to each independent component. These unmixing vectors form the unmixing matrix.

[0042] Step C5: Separate the whitened observation matrix into the unmixing matrix. Each of the statistically independent components contains [number] mutually independent components, and each independent component contains [number] components. Each time sampling point.

[0043] Step C6: Identify each independent component as an atmospheric delay independent component, a system noise independent component, and a deformation signal independent component based on their spatiotemporal characteristics. For the first... Calculate the following three quantitative indicators for each independent component: The first indicator is the spatial autocorrelation index. The Moran's I method is used to calculate and measure the spatial distribution regularity of this independent component. A binary spatial weight matrix is ​​used, with a weight of 1 when two pixels are adjacent and 0 otherwise. The value ranges from -1 to 1. The closer the value is to 1, the stronger the spatial positive autocorrelation, that is, the smoother the spatial distribution.

[0044] The second indicator is the spectral concentration. This measures the spectral energy distribution characteristics of the independent component. A discrete Fourier transform is performed on the independent component to calculate the normalized power spectral density. The ratio of the maximum value of the normalized power spectral density to its information entropy is taken as the spectral concentration. A larger value indicates that the spectral energy is more concentrated in a few frequency bands. The smaller the value, the more uniform the energy distribution of the spectrum.

[0045] The third indicator is the coefficient of determination for fitting the time trend. This measures the strength of the linear trend of the independent component over time. The least squares method is used to perform linear regression fitting on the time series of the independent components, and the ratio of the residual sum of squares to the total sum of squares is calculated to obtain the coefficients of determination. The value ranges from 0 to 1, and the larger the value, the more obvious the linear trend over time.

[0046] Based on the above three indicators, the following judgments are made in sequence, and each independent component is identified according to the first rule matched in the judgment order: The first step is to identify the independent components of the deformation signal: if and If it is identified as an independent component of the deformation signal, that is, the component has a clear linear trend in the time dimension and exhibits regional continuity in the spatial dimension, which is consistent with the physical characteristics of the surface deformation signal. The second step is to identify the atmospheric delayed independent components for the remaining independent components not identified in the first step: if and If it is identified as an atmospheric delay independent component, it means that the component has moderate or greater autocorrelation in the spatial dimension but lacks a significant temporal linear trend, which is consistent with the physical characteristics of atmospheric delay signals, which are large-scale smooth distribution in space and fluctuate with meteorological conditions in time. The third step involves uniformly identifying all remaining independent components not identified in the first and second steps as independent components of system noise, i.e., the spatial autocorrelation index of these components. And the coefficient of determination for fitting the time trend It tends to be randomly distributed in the spatial dimension and has no obvious trend in the temporal dimension, which is consistent with the statistical characteristics of system noise.

[0047] It should be noted that the above-mentioned sequential discrimination strategy ensures that each independent component has one and only one label result, the three categories are mutually exclusive and exhaustive, and cover all possible index value ranges.

[0048] In an embodiment of the present invention, the calculation results and discrimination process of the quantitative indicators of the five independent components are as follows: the fifth independent component , , The first step condition is met. and ), identified as independent components of the deformation signal; the first independent component , , The second step condition is met. and ), identified as the atmospheric delayed independent component; the third independent component , , The second step condition is met. and ), identified as the atmospheric delayed independent component; the second independent component , , If the conditions in steps one and two are not met, it is identified as an independent component of system noise by step three; the fourth independent component... , , If the conditions in the first and second steps are not met, it will be identified as an independent component of system noise by the third step.

[0049] Step 203, Denoising deformation phase sequence reconstruction; Independent components of the deformation signal are extracted and reconstructed to generate a denoised deformation phase sequence. The independent components identified as deformation signals are then subjected to an inverse transform using the corresponding column vectors in the mixing matrix, mapping them back to the original phase space to obtain the denoised deformation phase sequence. This sequence contains only the deformation signal component; atmospheric delay errors and system noise have been removed.

[0050] In an embodiment of the invention, after variational mode decomposition and independent component analysis, weak deformation signals in a large, deep landslide region were successfully separated. The average annual deformation rate of this region is approximately 3 to 8 millimeters, far below the detection threshold of traditional methods.

[0051] Step 204: Generation of multidimensional spatiotemporal feature tensors; Spatial registration is performed between the denoised deformation phase sequence and the set of terrain factors to generate a multidimensional spatiotemporal feature tensor, which includes the following steps: Step D1: Resample the denoised deformation phase sequence and terrain factor set to the same spatial resolution. Since synthetic aperture radar imagery and digital elevation model (DEM) may have different spatial resolutions, unification is necessary. In this embodiment, the denoised deformation phase sequence is resampled using bilinear interpolation, with the DEM resolution of 12.5 meters as the baseline. After resampling, the pixel positions of the two datasets correspond one-to-one.

[0052] Step D2: Expand the resampled denoised deformation phase sequence along the time dimension to obtain the deformation feature matrix. Assume the denoised deformation phase sequence contains... Each time sampling point, the target area includes If there are pixels, then the dimension of the deformation feature matrix is ​​. That is, each row corresponds to one pixel, and each column corresponds to the deformation phase value at a given time point. In an embodiment of the present invention, , , .

[0053] Step D3: Stack the terrain factor data from the terrain factor set along the feature channel dimension to obtain the terrain feature matrix. Assume the terrain factor set contains... If there are various terrain factors, then the dimension of the terrain feature matrix is... In this embodiment, It corresponds to six terrain factors: slope, aspect, plane curvature, profile curvature, topographic relief, and topographic roughness.

[0054] Step D4: Concatenate the deformation feature matrix and the terrain feature matrix along the feature channel dimension to generate a multidimensional spatiotemporal feature tensor. In this embodiment of the invention, the dimension of the multidimensional spatiotemporal feature tensor is... That is, 16 million pixels, each pixel has 221-dimensional features (215-dimensional deformation features and 6-dimensional terrain features).

[0055] Step 205, Construction of terrain topology map; A terrain topology map is constructed based on the spatial raster division of the target area. Multidimensional spatiotemporal feature tensors are then mapped to nodes in the terrain topology map to generate an attributed terrain feature map. This process includes the following steps: Step E1: Divide the target region into regular spatial grid cells. To reduce computational complexity while maintaining spatial resolution, a supercell method is used to aggregate adjacent pixels into grid cells. In this embodiment, every 4×4 adjacent pixels are aggregated into one grid cell, thus the target region is divided into... There are 16 raster cells. The feature vector of each raster cell is the average of the feature vectors of the 16 cells it contains.

[0056] Step E2: Define the set of nodes for the terrain topology map. Let the total number of grid cells be... Then the set of nodes Each node corresponds to a spatial grid cell. In an embodiment of the present invention, .

[0057] Step E3: Define the edge set of the terrain topology graph. Adjacent grid cells are connected by edges based on the eight-neighbor rule. The eight-neighbor rule is defined as follows: for a given grid cell located at the 8th ... Line number A column of grid cells has a neighborhood consisting of adjacent grid cells in eight directions: top, bottom, left, right, top-left, top-right, bottom-left, and bottom-right. Grid cells at boundaries are connected based on the actual number of their neighbors. If a node... With nodes If the corresponding grid cells are adjacent, an edge connection is established between them.

[0058] Step E4: Assign the feature vector corresponding to each grid cell in the multidimensional spatiotemporal feature tensor to the corresponding node of the terrain topology map. The initial feature vector of each node is the arithmetic mean of the feature vectors of all pixels within its corresponding grid cell. In this embodiment, each grid cell contains 16 pixels, therefore the initial feature vector of each node is the mean of the feature vectors of the 16 pixels, with a dimension of... A terrain feature map with attributes is constructed from the set of nodes, the set of edges, and the node feature matrix. The dimension of the node feature matrix is... .

[0059] Step 206, Graph attention network model processing; The attributed terrain feature map is input into a pre-trained graph attention network model, which outputs the hazard probability value and boundary confidence of each node.

[0060] The graph attention network model consists of an input layer, multiple graph attention convolutional layers, and an output layer. In this embodiment, the model contains three graph attention convolutional layers.

[0061] The input layer receives node feature vectors from an attributed terrain feature map and performs a linear transformation to reduce the feature dimension from... Dimensions are mapped to hidden dimensions.

[0062] Each graph attention convolutional layer processes node features through a graph attention mechanism. The specific process is as follows: First, linear transformations are performed on the features of the current node and its neighboring nodes. Then, the transformed features of the current node and its neighboring nodes are concatenated, and an attention score is calculated using the attention parameter vector activated by LeakyReLU. Next, softmax normalization is applied to the attention score to obtain attention weight coefficients. Finally, weighted aggregation is performed on the features of the neighboring nodes based on the attention weight coefficients, and the hidden state representation of the current node is updated using the ELU non-linear activation function. To enhance the model's expressive power, a multi-head attention mechanism is adopted. In this embodiment, there are 8 attention heads, and the outputs of each attention head are concatenated after independent calculation as the final output of the layer.

[0063] The output layer contains two parallel fully connected branches, which output the probability value of potential hazards for each node. and boundary confidence Both branches map the output to the 0-1 range using the Sigmoid function, where... Represents a node The probability of being a potential geological hazard point. Represents a node Confidence level at the boundary of the potential hazard.

[0064] The graph attention network model uses historical data containing labeled geological hazard hazard points for supervised learning during the training phase. Loss function. The weighted combination of binary cross-entropy loss and boundary-aware loss is used, and the calculation formula is as follows: ; in, For binary cross-entropy loss, For boundary-aware loss, These are the weighting coefficients.

[0065] Step 207: Determine the spatial location of potential geological hazard sites; Nodes with a hazard probability value greater than a preset probability threshold are selected to determine the spatial location of geological hazard hazard points. This process includes the following steps: Step F1: Traverse the hazard probability values ​​of each node in the terrain topology map. .

[0066] Step F2: Filter out potential hazards with a probability value greater than a preset probability threshold. Nodes that do not meet the criteria are marked as potential hazard nodes. The preset probability threshold is pre-calculated and can be adjusted according to actual needs (such as considering both false negative and false positive rates). In this embodiment, Nodes with a probability value greater than 0.7 are marked as potential hazard nodes.

[0067] Let the set of candidate nodes for potential hazards be . In an embodiment of the present invention, 857 potential hazard candidate nodes were obtained through screening.

[0068] Step F3: Obtain the center coordinates of the spatial grid cell corresponding to the candidate node of the geological hazard to determine the spatial location of the geological hazard point. Let the grid cell be located at the... Line number Column, grid cell side length is The coordinates of the top left corner of the target area are Then the center coordinates of the grid cell for: ; ; The pixel coordinates are converted to a geographic coordinate system to output the spatial location of potential geological hazard points. In this embodiment, the WGS-84 geographic coordinate system is used, and the output format is latitude and longitude coordinates.

[0069] Step 208, Region Growth Boundary Tracing; Based on the boundary confidence level, the region growth boundary tracking is performed on adjacent hazard nodes to generate hazard body boundary vector data.

[0070] Step G1: Select the candidate node with the highest probability value of the hidden danger as the seed node.

[0071] Step G2: Initialize the current potential hazard area, initialize the boundary node queue, and add the neighboring nodes of the seed node to the boundary node queue.

[0072] Step G3: Retrieve neighboring nodes from the boundary node queue The boundary confidence level is used to determine whether a neighboring node belongs to the same potential hazard. The determination criterion is the boundary confidence level of the neighboring node. Greater than the preset boundary confidence threshold Furthermore, the neighboring node belongs to the set of potential hazard candidate nodes; in this embodiment .

[0073] If the conditions are met, then the neighboring nodes will be... Included in the current potential hazard area, and Unvisited neighboring nodes are added to the boundary node queue.

[0074] Step G4: Iterate through step G3 until the boundary node queue is empty, meaning the boundary of the current potential hazard area has converged.

[0075] Step G5: Extract the outer contour of the current hazard area. Using an eight-neighbor boundary tracing algorithm, starting from the top-left corner node of the hazard area, trace the boundary nodes clockwise until returning to the starting node. The sequence of boundary nodes constitutes the closed boundary of the hazard.

[0076] Step G6: Convert the coordinate sequence of the boundary nodes into a vector polygon format to generate the boundary vector data of the potential hazard. The vector data format adopts the GeoJSON standard and includes geometric type, coordinate array, and attribute information.

[0077] Step G7: Remove the processed hazard candidate nodes from the hazard candidate node set. If the set is not empty, return to step G1 to process the next hazard.

[0078] In the embodiments of the present invention, a total of three independent geological hazard bodies were identified, with areas of approximately 0.8 square kilometers, 1.2 square kilometers, and 0.5 square kilometers, respectively.

[0079] Remote access terminal 3 is configured to receive spatial location and boundary vector data of geological hazard points and hazard bodies transmitted by data processing center 2, and perform visualization and storage processing. The specific execution process includes the following steps: Step 301, generating a visualization layer for the distribution of potential hazards, specifically includes the following steps: Step H1: Load the preset Geographic Information System (GIS) base map. The GIS base map includes a topographic base map and an administrative division layer. The topographic base map is composed of a mountain shadow map and a contour map generated by digital elevation model rendering, overlaid together. The administrative division layer contains the administrative boundaries at the provincial, municipal, county, and township levels.

[0080] Step H2: Convert the spatial locations of geological hazard hazard points into a point feature layer. Each hazard point is represented as a point feature, and its attributes include: point number, longitude coordinates, latitude coordinates, hazard probability value, and identification time. Point features are displayed using red circular symbols, and the size of the symbols is proportional to the hazard probability value.

[0081] Step H3: Convert the boundary vector data of the hazard bodies into a polygon feature layer. Each hazard body is represented as a polygon feature, and the feature attributes include: hazard body number, area, perimeter, average hazard probability, and average boundary confidence. Polygon features are displayed with a semi-transparent red fill and a dark red border.

[0082] Step H4: Overlay the point feature layer and the polygon feature layer onto the geographic information system base map according to spatial coordinates. The layer overlay order from bottom to top is: topographic base map, administrative division layer, polygon feature layer, and point feature layer. The generated hazard distribution visualization layer supports interactive operations such as zooming, panning, layer control, and feature querying.

[0083] Step 302, Visualize layer storage; The visualization layer showing the distribution of potential hazards is stored in a spatial database. The spatial database uses a relational database that supports spatial indexing. The data table structure includes: Hazard Point Table: Fields include point number (primary key), geometric object (point type), hazard probability value, identification time, and area code; Hazard body table: Fields include hazard body number (primary key), geometric object (face type), area, perimeter, average hazard probability, average boundary confidence, identification time, and region code; Identification record table: Fields include record number (primary key), identification time, processing parameters, number of hidden danger points, and number of hidden danger bodies.

[0084] Spatial indexes are created to accelerate spatial queries. The spatial indexes use an R-tree index structure and support range queries, proximity queries, and spatial relationship queries.

[0085] In embodiments of the present invention, after the identification results are stored in a spatial database, they can be queried and visualized via a remote access terminal 3. Disaster prevention and mitigation management personnel can formulate monitoring and early warning plans and engineering remediation measures based on the distribution of potential hazards.

[0086] The specific execution flow of this invention is as follows: Step 100: Data acquisition terminal 1 acquires a multi-temporal synthetic aperture radar interferometric image sequence of the target area, performs spatiotemporal baseline registration on the multi-temporal synthetic aperture radar interferometric image sequence, and generates a temporal interferometric phase dataset; acquires digital elevation model data of the target area, and calculates a set of terrain factors based on the digital elevation model data; Step 200: Data Processing Center 2 receives the temporal interferometric phase dataset and the topographic factor set. It performs variational mode decomposition on the temporal interferometric phase dataset pixel by pixel, decomposing the temporal interferometric phase of each pixel into multiple intrinsic mode components. An observation matrix is ​​constructed for these intrinsic mode components, and independent component analysis is performed on the observation matrix to separate the atmospheric delay independent component, system noise independent component, and deformation signal independent component. The deformation signal independent component is extracted and reconstructed to generate a denoised deformation phase sequence. Spatial registration is performed between the denoised deformation phase sequence and the topographic factor set to generate a multidimensional spatiotemporal feature tensor. A topographic topology map is constructed based on the spatial grid division of the target area, and the multidimensional spatiotemporal feature tensor is mapped to the nodes of the topographic topology map to generate an attributed topographic feature map. The attributed topographic feature map is input into a pre-trained graph attention network model, which outputs the hazard probability value and boundary confidence of each node. Nodes with hazard probability values ​​greater than a preset probability threshold are selected to determine the spatial location of geological hazard points. Regional growth boundary tracking is performed on adjacent hazard nodes based on the boundary confidence to generate hazard body boundary vector data. Step 300: Remote access terminal 3 receives spatial location data of geological hazard points and boundary vector data of hazard bodies; overlays the spatial location data of geological hazard points and boundary vector data of hazard bodies onto the geographic information system base map to generate a hazard distribution visualization layer; and stores the hazard distribution visualization layer into the spatial database.

[0087] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A device for precise location of geological hazard hazard points based on intelligent remote sensing interpretation combined with topographic factors, characterized in that: include: Data acquisition terminals, data processing centers, and remote access terminals; The data acquisition terminal is used to perform the following steps: Acquire multi-temporal synthetic aperture radar (SAR) interferometric image sequences of the target area, perform spatiotemporal baseline registration on the multi-temporal SAR SAR interferometric image sequences, and generate a temporal interferometric phase dataset; acquire digital elevation model (DEM) data of the target area, and calculate a set of terrain factors based on the DEM data; The data processing center is used to perform the following steps: Receive a temporal interferometric phase dataset and a set of terrain factors, and perform variational mode decomposition on the temporal interferometric phase dataset pixel by pixel, decomposing the temporal interferometric phase of each pixel into multiple intrinsic mode components; An observation matrix is ​​constructed for multiple intrinsic mode components. Independent component analysis is performed on the observation matrix to separate the atmospheric delay independent component, the system noise independent component, and the deformation signal independent component. The deformation signal independent component is extracted and reconstructed to generate a denoised deformation phase sequence. Spatial registration is performed between the denoised deformation phase sequence and the set of topographic factors to generate a multidimensional spatiotemporal feature tensor. A topographic topology map is constructed based on the spatial grid division of the target area. Multidimensional spatiotemporal feature tensors are mapped to the nodes of the topographic topology map to generate attributed topographic feature maps. The attributed topographic feature maps are input into a pre-trained graph attention network model, which outputs the hazard probability value and boundary confidence of each node. Nodes with hazard probability values ​​greater than a preset probability threshold are selected to determine the spatial location of geological hazard points. Based on the boundary confidence, regional growth boundary tracking is performed on adjacent hazard nodes to generate hazard body boundary vector data. The remote access terminal is used to perform the following steps: Receive spatial location data of geological hazard hazard points and vector data of hazard body boundaries; overlay the spatial location data of geological hazard hazard points and vector data of hazard body boundaries onto the geographic information system base map to generate a hazard distribution visualization layer; store the hazard distribution visualization layer into a spatial database.

2. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors as described in claim 1, characterized in that, The generation of the temporal interferometric phase dataset includes: selecting the master image from the multi-temporal synthetic aperture radar interferometric image sequence and using the remaining images as auxiliary images; calculating the offset between images using the intensity cross-correlation method, establishing the geometric mapping relationship between the auxiliary images and the master image through polynomial fitting, resampling the auxiliary images into the geometric framework of the master image to eliminate the geometric offset between images; calculating the temporal baseline and spatial baseline between each auxiliary image and the master image, where the temporal baseline is defined as the difference between the acquisition time of the auxiliary image and the acquisition time of the master image, and the spatial baseline is defined as the projection distance between the satellite orbit positions of the auxiliary image and the master image in the direction perpendicular to the radar line of sight; selecting image pairs that meet the conditions of the spatial baseline being less than two-thirds of a preset critical baseline and the temporal baseline being less than a preset time decoherence threshold to construct an interferometric pair set; multiplying the complex values ​​of the master image and the auxiliary image by their conjugates to obtain the complex interferometric map, extracting the phase information from the complex interferometric map to obtain the entangled interferometric phase, and using differential interferometry to sequentially remove the flat phase and terrain phase from the original interferometric phase to obtain the differential interferometric phase map; arranging each interferometric phase map in chronological order to generate the temporal interferometric phase dataset.

3. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 2, characterized in that, The method for calculating a set of terrain factors based on digital elevation model (DEM) data includes: performing a sliding window traversal on the DEM data to calculate the slope and aspect data for each pixel, where the slope data represents the rate of change of elevation and the aspect data represents the azimuth of the direction of maximum slope; calculating the planar curvature and profile curvature data for each pixel based on the second-order partial derivatives of the DEM data; calculating the difference between the maximum and minimum elevations within a preset neighborhood window to obtain terrain relief data; calculating the standard deviation of elevation values ​​within the neighborhood window to obtain terrain roughness data; and combining the slope data, aspect data, planar curvature data, profile curvature data, terrain relief data, and terrain roughness data into a set of terrain factors.

4. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 3, characterized in that, Performing variational mode decomposition on a pixel-by-pixel basis on a temporal interferometric phase dataset includes: constructing a constrained optimization problem for variational mode decomposition, taking the phase time series of each pixel in the temporal interferometric phase dataset as the signal to be decomposed; iteratively solving the constrained optimization problem using the alternating direction multiplier method to obtain the center frequency and bandwidth of each intrinsic mode component; and performing adaptive decomposition on the phase time series based on the center frequency and bandwidth to output K intrinsic mode components, where K is a positive integer.

5. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 4, characterized in that, Performing independent component analysis on the observation matrix includes: constructing the observation matrix by arranging the K intrinsic mode components in rows; performing centering and whitening processes on the observation matrix in sequence to obtain a whitened matrix; performing independent component analysis on the whitened matrix based on the non-Gaussianity maximization criterion to obtain an unmixed matrix; separating the observation matrix into statistically independent components through the unmixed matrix; and identifying each independent component as an atmospheric delay independent component, a system noise independent component, and a deformation signal independent component based on the spatiotemporal characteristics of each component.

6. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 5, characterized in that, The generation of the multidimensional spatiotemporal feature tensor includes: resampling the denoised deformation phase sequence and the set of terrain factors to the same spatial resolution; expanding the resampled denoised deformation phase sequence along the time dimension to obtain the deformation feature matrix; stacking the terrain factor data of each terrain factor in the set along the feature channel dimension to obtain the terrain feature matrix; and concatenating the deformation feature matrix and the terrain feature matrix along the feature channel dimension to generate the multidimensional spatiotemporal feature tensor.

7. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 6, characterized in that, The generation of attributed terrain feature maps includes: dividing the target area into regular spatial grid cells; defining a set of nodes for the terrain topology map, with each node corresponding to a spatial grid cell; defining a set of edges for the terrain topology map, with adjacent grid cells connected by an eight-neighborhood connection rule; and assigning the feature vector corresponding to each grid cell in the multidimensional spatiotemporal feature tensor to the corresponding node of the terrain topology map to generate an attributed terrain feature map.

8. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 7, characterized in that, The attributed terrain feature map is input into a pre-trained graph attention network model, which outputs the hazard probability value and boundary confidence of each node. The graph attention network model includes an input layer, multiple graph attention convolutional layers, and an output layer. The input layer receives the node feature vectors of the attributed terrain feature map. Each graph attention convolutional layer calculates the attention weight coefficients between the current node and its neighboring nodes through a multi-head attention mechanism, performs weighted aggregation on the features of neighboring nodes based on the attention weight coefficients, and updates the hidden state representation of the current node. The output layer performs classification mapping on the hidden state representation of the nodes output by the last graph attention convolutional layer through two parallel fully connected branches, and outputs the hazard probability value and boundary confidence of each node.

9. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 8, characterized in that, Based on the boundary confidence score, the region growth boundary tracking is performed on adjacent hazard nodes to generate hazard body boundary vector data. This includes: selecting the hazard candidate node with the highest hazard probability value as the seed node; starting from the seed node, determining whether neighboring nodes belong to the same hazard body based on the boundary confidence score, and including neighboring nodes that meet the boundary confidence score condition into the current hazard body region; iteratively performing region growth until the boundary of the current hazard body region converges; extracting the boundary node sequence of the outer contour of the current hazard body region, converting the coordinate sequence of the boundary nodes into vector polygon format, and generating hazard body boundary vector data.

10. The device for precise location of geological hazard points based on intelligent remote sensing interpretation combined with topographic factors according to claim 9, characterized in that, The criteria for determining whether a neighboring node belongs to the same hazard body based on boundary confidence are as follows: if the boundary confidence of a neighboring node is greater than the preset boundary confidence threshold, and the neighboring node belongs to the set of hazard candidate nodes, then the neighboring node belongs to the same hazard body; otherwise, the neighboring node does not belong to the same hazard body.

Images