Slope instability early warning method and system based on unmanned aerial vehicle remote sensing data

Abstract

The application discloses a slope instability early warning method and system based on unmanned aerial vehicle remote sensing data, and relates to the technical field of slope instability early warning. In view of the problem that the slope instability precursor recognition accuracy is not high, leading to lag of early warning and even complete failure, the application comprises the following steps: acquiring three-dimensional point cloud data of a target slope, constructing a dynamic space-time graph network, constructing structure evolution indexes of each point cloud node, and determining a data fidelity stiffness parameter in a preset local energy functional equation; updating the local energy functional equation based on the data fidelity stiffness parameter, solving the minimum value of the updated local energy functional equation, and outputting real ground elevation data covering a smoothing processing area and an exempted reserved area; reconstructing a high-fidelity digital elevation model, and performing data space difference calculation on high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data, and performing slope instability early warning. The application improves the accuracy and robustness of slope instability precursor recognition.

Description

Technical Field

[0001] This invention relates to the field of slope instability early warning technology, specifically to a slope instability early warning method and system based on UAV remote sensing data. Background Technology

[0002] In actual natural mountain or complex surface environments, the early signs of slope instability often manifest as tiny cracks or bulges on the surface, with extremely small spatial scales and high spatial frequencies. This severe overlap of surface physical precursor signals with environmental noise at spatial frequencies constitutes a technical bottleneck that is difficult to overcome with current technology.

[0003] Because existing point cloud filtering algorithms generally rely on single spatial scale features or globally fixed geometric smoothing parameters for noise reduction, they cannot effectively decouple effective precursor deformations and environmental noise in the same high-frequency band within a single static spatial domain when faced with the aforementioned frequency band overlap phenomenon. This leads to existing algorithms, when performing high-frequency noise filtering, inevitably treating real, minute precursor deformations as spatial outliers for forced smoothing in mathematical processing. This traditional static filtering mechanism produces a severe "oversampling" phenomenon in complex terrain, indiscriminately smoothing out key micro-fracture features in the early stages of slope evolution, causing serious distortion in subsequent deformation calculations, ultimately resulting in delayed or even complete failure of early warning for slope instability.

[0004] In view of the above, this application is hereby submitted. Summary of the Invention

[0005] The technical problem this invention aims to solve is that existing slope instability early warning methods fail to distinguish between micro-fracture characteristics and environmental disturbance noise (such as vegetation noise), and do not achieve high-sensitivity characterization of the evolution process of slope micro-fractures. This results in low accuracy in identifying early signs of slope instability, leading to delayed or even completely ineffective early warnings. The purpose of this invention is to provide a slope instability early warning method and system based on UAV remote sensing data. This method achieves high-sensitivity characterization of the evolution process of slope micro-fractures, improving the accuracy and robustness of identifying early signs of slope instability. Furthermore, it enables high-precision graded early warning of slope instability risk, solving the problems of delayed or completely ineffective early warnings.

[0006] This invention is achieved through the following technical solution:

[0007] In a first aspect, the present invention provides a slope instability early warning method based on UAV remote sensing data, the method comprising:

[0008] Acquire spatially registered 3D point cloud data of the target slope at a preset number of continuous time nodes. The 3D point cloud data includes the local normal vector of each point cloud node and the basic feature vector of the original elevation data.

[0009] Based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data, a dynamic spatiotemporal graph network is constructed.

[0010] Based on the dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vector to obtain the structural evolution index of each point cloud node.

[0011] Based on the structural evolution index, determine the data fidelity stiffness parameters in the preset local energy functional equation;

[0012] The local energy functional equation is updated based on the data-fidelity stiffness parameter. The minimum value of the updated local energy functional equation is solved by a preset solution algorithm, and the output is real surface elevation data covering the smoothed area and the exempted retention area.

[0013] Based on multiple sets of real surface elevation data, high-fidelity digital elevation models corresponding to each time node are reconstructed, and data spatial difference calculations are performed on the high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data.

[0014] In pure three-dimensional displacement field data, if the target deformation rate in a region where the structural evolution index is greater than the preset evolution threshold exceeds the preset instability critical value, slope instability early warning data will be output to the terminal.

[0015] Furthermore, the steps for obtaining the local normal vectors of each point cloud node and the fundamental feature vectors of the original elevation data are as follows:

[0016] At the current time point, extract the three-dimensional spatial coordinates of the point cloud nodes and generate neighboring point set data based on the preset neighborhood radius;

[0017] A local covariance matrix is ​​constructed based on the three-dimensional spatial coordinates of adjacent point set data. Eigenvalue decomposition is performed on the local covariance matrix, and the eigenvectors are extracted as local normal vectors.

[0018] Extract the difference in three-dimensional spatial coordinates of point cloud nodes between the current time node and the previous time node, and divide the difference in three-dimensional spatial coordinates by the corresponding time interval data to obtain the instantaneous displacement rate;

[0019] The basic feature vector is generated by splicing together the original elevation data corresponding to the three-dimensional spatial coordinates, the local normal vector data, and the instantaneous displacement rate.

[0020] Furthermore, based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data, a dynamic spatiotemporal graph network is constructed, including:

[0021] Get any two point cloud nodes at the same time point and calculate the inner product of the two point cloud nodes.

[0022] Calculate the spatial edge weights based on the inner product and the Euclidean distance between any two point cloud nodes.

[0023] Connect point cloud nodes representing the same spatial physical location at adjacent time points to generate time edge weights;

[0024] Point cloud nodes are used as node feature inputs, and a global spatiotemporal graph matrix is ​​established by combining spatial edge weights and temporal edge weights to generate a dynamic spatiotemporal graph network.

[0025] Furthermore, the specific steps for generating a dynamic spatiotemporal graph network are as follows:

[0026] Extract a preset number of continuous time nodes ( Given the total number of consecutive time nodes, all point cloud nodes participating in the construction of the graph network are represented by N, with a total number of point cloud nodes in a single period. The basic feature vector C corresponding to each point cloud node is arranged according to its spatial index and time sequence, generating a dimension of... The multidimensional node feature matrix X; where F is the length of the basic feature vector C;

[0027] Traverse all node pairs at any given time point and construct a spatial adjacency submatrix based on the spatial edge weights of the node pairs. The spatial adjacency submatrix is ​​used to represent the spatial topological connectivity of the slope surface at each time point in the graph network.

[0028] Based on the connection relationship between any point cloud node at the same spatial physical location at time node t and time node t+1, construct a time adjacency submatrix;

[0029] The spatial adjacency submatrix and the temporal adjacency submatrix are concatenated in blocks to construct a global spatiotemporal graph matrix. ;

[0030] The multidimensional node feature matrix X and the global spatiotemporal graph matrix Binding and encapsulation are performed to generate a dynamic spatiotemporal graph network. .

[0031] Furthermore, based on the dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vectors to obtain the structural evolution indicators of each point cloud node, including:

[0032] By using the spatial convolutional layer in the dynamic spatiotemporal graph network and aggregating the neighborhood node features of the basic feature vector according to the spatial edge weights, the spatial anisotropic feature data of the generated point cloud nodes are extracted.

[0033] By using the time-series convolutional layer in the dynamic spatiotemporal graph network and performing time-series gradient integration on the instantaneous displacement rate in the basic feature vector based on the time edge weights, the temporal monotonic cumulative feature data of the generated point cloud nodes is extracted.

[0034] Normalized multiplication is performed on spatially anisotropic feature data and temporally monotonic cumulative feature data to output a structural evolution index with numerical constraints within a preset numerical range.

[0035] Furthermore, based on structural evolution indices, the data fidelity stiffness parameters in the preset local energy functional equations are determined, including:

[0036] If the structural evolution index is greater than the preset evolution threshold, the data fidelity stiffness parameter will be mapped to the preset maximum value and the smoothing exemption mechanism for micro-fracture data will be triggered; otherwise, the current structural evolution index will be mapped to the data fidelity stiffness parameter.

[0037] Mapping data fidelity stiffness parameters to preset maxima includes:

[0038] Read the preset maximum value stored in the preset database. The magnitude of the preset maximum value is greater than the preset smoothing weight corresponding to the spatial smoothing term in the local energy functional equation.

[0039] If the structural evolution index is greater than the preset evolution threshold, a data assignment operation is performed, and the preset maximum value is assigned to the data fidelity stiffness parameter.

[0040] The assigned data fidelity stiffness parameter is substituted into the data fidelity term of the local energy functional equation, and in each iteration of the preset solution algorithm, the elevation residual data between the real surface elevation data to be solved and the original elevation data is calculated.

[0041] Numerical penalty calculations are applied to the elevation residual data based on the assigned data fidelity stiffness parameters, and when the energy value of the local energy functional equation converges to a minimum, the output value of the true surface elevation data is directly constrained to be equal to the original elevation data.

[0042] Furthermore, after mapping the data fidelity stiffness parameter to a preset maximum value, it also includes:

[0043] If the structural evolution index is less than the preset denoising threshold, and the preset denoising threshold is less than the preset evolution threshold, then the data fidelity stiffness parameter is mapped to the preset minimum value, so that the real surface elevation data is close to the macroscopic low-frequency elevation data output by the spatial smoothing term, triggering a deep denoising mechanism for vegetation cover data.

[0044] If the structural evolution index is between the preset denoising threshold and the preset evolution threshold, then linear interpolation is performed based on the current value of the structural evolution index, and the output is a data-fidelity stiffness parameter with a linear transition.

[0045] Furthermore, the minimum value of the updated local energy functional equation is solved using a preset solution algorithm, including:

[0046] Substitute the original elevation data and data fidelity stiffness parameters into the data fidelity term to establish the first polynomial describing the data fidelity error.

[0047] Substitute the actual surface elevation data to be solved and the spatial neighborhood elevation data of the point cloud nodes into the spatial smoothing term to establish a second polynomial describing the surface smoothing error.

[0048] Combine the first polynomial and the second polynomial to transform them into the objective functional matrix;

[0049] Based on a preset sparse matrix optimization algorithm, the target functional matrix is ​​differentiated and iteratively degraded. When the iteration residual is less than the preset convergence threshold, the real surface elevation data that satisfies the energy minimization condition is calculated.

[0050] Furthermore, spatial difference calculations are performed on the high-fidelity digital elevation models at adjacent time points to generate pure three-dimensional displacement field data, including:

[0051] Extract the first digital elevation model data corresponding to the first time node from the pre-stored high-fidelity digital elevation model, and extract the second digital elevation model data corresponding to the second time node adjacent to the first time node;

[0052] Based on the preset spatial grid resolution, the first digital elevation model data and the second digital elevation model data are subjected to spatial grid projection and alignment processing to generate aligned grid data with the same plane coordinate sequence.

[0053] Iterate through each plane coordinate point in the aligned grid data, and extract the first elevation data of that plane coordinate point in the first digital elevation model data and the second elevation data in the second digital elevation model data.

[0054] The elevation change data for each plane coordinate point is calculated by subtracting the second elevation data from the first elevation data.

[0055] The calculated elevation change data is extracted, and the elevation change data is combined with the corresponding plane coordinate points to perform a three-dimensional vector stitching operation to generate pure three-dimensional displacement field data containing spatial coordinate information and vertical deformation values.

[0056] Secondly, the present invention provides a slope instability early warning system based on UAV remote sensing data, the system comprising:

[0057] The data acquisition module is used to acquire three-dimensional point cloud data of the target slope after spatial registration at a preset number of continuous time nodes. The three-dimensional point cloud data includes the local normal vector of each point cloud node and the basic feature vector of the original elevation data.

[0058] The feature processing module is used to construct a dynamic spatiotemporal graph network based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data; based on the dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vectors to obtain the structural evolution index of each point cloud node.

[0059] The filtering control module is used to determine the data fidelity stiffness parameters in the preset local energy functional equation based on the structural evolution index; update the local energy functional equation based on the data fidelity stiffness parameters; solve the minimum value of the updated local energy functional equation through a preset solution algorithm; and output the real surface elevation data covering the smoothed area and the exempted retention area.

[0060] The early warning output module is used to reconstruct high-fidelity digital elevation models corresponding to each time node based on multiple sets of real surface elevation data, and to perform data spatial difference calculation on the high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data. In the pure three-dimensional displacement field data, if the target deformation rate in the area where the structural evolution index is greater than the preset evolution threshold exceeds the preset instability critical value, the slope instability early warning data is output to the terminal.

[0061] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0062] 1. This invention relates to a slope instability early warning method and system based on UAV remote sensing data. By constructing a dynamic spatiotemporal graph network and integrating the consistency of local normal vector angles and instantaneous displacement rate characteristics, this invention achieves a highly sensitive characterization of the evolution process of slope micro-fractures. Unlike traditional methods that rely on single geometric features or temporal differences, this invention constructs spatial edge weights based on the inner product of local normal vectors, which can more accurately capture the shear slip trend of internal structural surfaces of soil and rock masses. At the same time, combined with the time monotonic cumulative features generated by the gradient integral of the time series, it helps to amplify the nonlinear deformation signal of the micro-fracture region. The resulting structural evolution index can not only quantify the joint characteristics of point cloud nodes in terms of anisotropy and temporal monotonicity, but also provide a physically meaningful adaptive weight basis for subsequent energy functional filtering, which is beneficial to improving the accuracy and robustness of slope instability precursor identification.

[0063] 2. This invention relates to a slope instability early warning method and system based on UAV remote sensing data. The invention introduces a local energy functional filtering mechanism driven by structural evolution indices, achieving differentiated processing of the original elevation data and possessing the dual advantages of suppressing strong noise and preserving micro-fracture features. Furthermore, by mapping structural evolution indices to data-fidelity stiffness parameters and coordinating with a preset evolution threshold to trigger an exemption mechanism for smoothing, maximum stiffness constraints can be automatically applied to suspected fracture areas, ensuring that the actual surface elevation data closely approximates the original observations in numerical accuracy, thereby avoiding excessive smoothing of key deformation information. Simultaneously, for areas with strong noise such as vegetation cover, a preset denoising threshold triggers minimum stiffness mapping, smoothing the filtering results to the macroscopic topographic background and effectively filtering out non-target disturbances. Based on this adaptive stiffness adjustment mechanism, the filtered digital elevation model maintains the integrity of the original micro-topography while improving the purity and reliability of subsequent displacement field inversion.

[0064] 3. This invention relates to a slope instability early warning method and system based on UAV remote sensing data. By generating a pure three-dimensional displacement field and combining it with structural evolution indicators for regional constraint determination, this invention achieves high-precision graded early warning of slope instability risk. Simultaneously, spatial difference based on a high-fidelity digital elevation model between time nodes effectively eliminates systematic errors such as atmospheric disturbances and vegetation obstruction, outputting only a three-dimensional displacement vector reflecting actual terrain changes. At the early warning logic level, this invention uses regions where the deformation rate exceeds a preset critical value and the structural evolution indicator is greater than the evolution threshold as a joint criterion, avoiding false alarms or missed alarms caused by a single indicator and improving the field adaptability of the early warning system. Attached Figure Description

[0065] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0066] Figure 1 The present invention provides a flowchart for a slope instability early warning method based on UAV remote sensing data. Figure 1 ;

[0067] Figure 2 The present invention provides a flowchart for a slope instability early warning method based on UAV remote sensing data. Figure 2 ;

[0068] Figure 3 This is a structural block diagram of the slope instability early warning system based on UAV remote sensing data of the present invention. Detailed Implementation

[0069] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.

[0070] like Figure 1 As shown, the present invention includes: First, an unmanned aerial vehicle (UAV) conducts multiple aerial surveys of a target slope according to a preset flight plan, acquiring a batch of raw 3D point cloud data in each survey. Spatial registration of the point clouds from each phase is performed using the Iterative Closest Point (ICP) algorithm or Global Navigation Satellite System (GNSS) control points to eliminate UAV flight attitude errors and coordinate system deviations, resulting in target 3D point cloud data with unified spatial coordinates. Then, for each point cloud node in the target 3D point cloud data, local geometric features are constructed in its neighborhood, local normal vectors and instantaneous displacement rates are calculated, and these are concatenated to generate a basic feature vector, which serves as input for subsequent graph network inference. Next, spatial edges are constructed based on the directional consistency of local normal vectors between adjacent point cloud nodes, and temporal edges are constructed based on the temporal correspondence of the same physical location, forming a dynamic spatiotemporal graph network. The basic feature vector is input into this network, and structural evolution indicators for each point cloud node are extracted through graph convolution operations. Then, using structural evolution indices as the core, the data-fidelity stiffness parameters in the local energy functional equations are adaptively adjusted. For nodes where the structural evolution indices exceed a preset evolution threshold, a smoothing exemption mechanism is triggered, forcibly preserving the original elevation data. For the remaining nodes, the stiffness parameters are smoothly mapped using a tangent function model to achieve adaptive filtering. Subsequently, using the real surface elevation data from each period output after adaptive filtering, a high-fidelity digital elevation model is reconstructed, and differential calculations are performed on the models from adjacent periods to obtain pure three-dimensional displacement field data. Finally, in the pure three-dimensional displacement field data, the target deformation rate in the region where the structural evolution indices exceed the threshold is judged. If it exceeds a preset instability critical value, slope instability early warning data is output to the terminal.

[0071] Example 1

[0072] like Figure 2 As shown, the present invention provides a slope instability early warning method based on UAV remote sensing data, the method comprising:

[0073] S1, acquire the three-dimensional point cloud data of the target slope after spatial registration at a preset number of continuous time nodes. The three-dimensional point cloud data includes the local normal vector of each point cloud node and the basic feature vector C of the original elevation data.

[0074] In this embodiment, the original three-dimensional point cloud data refers to the set of point coordinates of the ground surface and its cover (vegetation, gravel, etc.) in three-dimensional space, which is obtained by a UAV equipped with a LiDAR sensor or a photogrammetry system (such as Structure from Motion, SfM). Each point has three-dimensional spatial coordinates, and some schemes also include attribute information such as reflection intensity.

[0075] In this context, a point cloud node refers to each discrete spatial sampling point with three-dimensional coordinates in the three-dimensional point cloud data, which is the basic processing unit for graph network construction and feature extraction in this invention.

[0076] Among them, the local normal vector refers to the unit normal vector describing the orientation of the local surface where the point cloud node is located. It is obtained by performing principal component analysis (PCA) on the point set in the neighborhood of the node and is used to characterize the geometric features of the local surface.

[0077] Step S1 above is the data preparation and feature representation stage of the present invention. Since the evolution of slope instability is a gradual physical process spanning time, static point cloud data at a single time node can only reflect the spatial morphology of the slope at that specific time point and cannot reveal its evolutionary trend over time. Therefore, it is necessary to acquire multi-period point cloud data at continuous time nodes and eliminate measurement system errors through spatial registration to ensure the comparability of data from different periods. Furthermore, the subsequent construction and inference of the dynamic spatiotemporal graph network rely on the structured feature representation of each point cloud node. Since the original three-dimensional coordinates are insufficient to carry enough local geometric and temporal information, it is necessary to further extract basic feature vectors, including local normal vectors and instantaneous displacement rates, as the input basis for graph network inference.

[0078] Taking a mountain slope under monitoring as an example, the UAV conducts an aerial survey flight on the target slope every 7 days to continuously acquire 5 periods of original 3D point cloud data (i.e., the total number of consecutive time nodes is preset to 5). After spatial registration, each period of point cloud contains about 5 million point cloud nodes in a unified coordinate system. Each node carries 3D spatial coordinate information, providing a data foundation for the subsequent extraction of basic feature vectors.

[0079] In this embodiment, the steps for obtaining the local normal vector of each point cloud node and the basic feature vector of the original elevation data are as follows:

[0080] S11. At the current time point, extract the three-dimensional spatial coordinates of the point cloud nodes and generate neighboring point set data according to the preset neighborhood radius; wherein, the preset neighborhood radius is a preset value, which is preset by preset staff based on historical data and experience rules, and stored in a preset database.

[0081] S12, construct a local covariance matrix based on the three-dimensional spatial coordinates of adjacent point set data, perform eigenvalue decomposition on the local covariance matrix, and extract the eigenvectors as local normal vectors;

[0082] Specifically, local normal vectors are extracted using Principal Component Analysis (PCA): Centered on a point cloud node, neighboring point sets are searched within a preset neighborhood radius (typically 5 to 10 times the average point spacing in the point cloud; for example, a neighborhood radius of 0.5m is set when the point spacing is approximately 0.1m), and a local covariance matrix C is constructed. cov : Where T is the transpose operator and k is the number of points in the neighborhood. Let be the three-dimensional spatial coordinates of the i-th neighboring point. Let C be the mean coordinates of the neighborhood point set; for the local covariance matrix C cov Eigenvalue decomposition yields three eigenvalues. and the corresponding feature vectors Minimum eigenvalue Corresponding feature vector This refers to the data of the local normal vector, which physically represents the direction of the normal to the local surface formed by the set of neighborhood points. Each node should have at least 15 neighboring points in its neighborhood to ensure the numerical stability of the local covariance matrix. When the number of neighborhood points is insufficient, the neighborhood radius can be increased.

[0083] S13, extract the difference in three-dimensional spatial coordinates of point cloud nodes between the current time node and the previous time node, divide the difference in three-dimensional spatial coordinates by the corresponding time interval data to obtain the instantaneous displacement rate;

[0084] Specifically, instantaneous displacement rate The calculation formula is: ,in Let be the three-dimensional spatial coordinates of the point cloud node at the current time node t. Δt represents the three-dimensional spatial coordinates of the point cloud node at the previous time node t-1 (the correspondence is determined by nearest neighbor matching), and Δt is the time interval between two adjacent time nodes.

[0085] S14: Based on the original elevation data, local normal vector data, and instantaneous displacement rate corresponding to the three-dimensional spatial coordinates, a basic feature vector is generated by splicing them together.

[0086] Specifically, the original elevation data z and the three-dimensional data of the local normal vector are... and instantaneous displacement rate 3D vector By concatenating the components, we obtain the basic feature vector C of length F=7. For example, taking the slope monitoring scenario mentioned above, suppose the coordinates of a point cloud node P in period 4 (the current time node) are (100.000m, 200.000m, 450.031m), and the corresponding coordinates in period 3 (the previous time node) are (100.000m, 200.000m, 450.012m). The time interval between two adjacent time nodes is 7 days. Then, the instantaneous displacement rate in... This indicates that the point has a continuous upward trend in the vertical direction. The local normal vector [0.02, 0.03, 0.999] obtained by joint neighborhood PCA calculation (approximately vertically upward, indicating that the local surface is relatively flat) is finally generated as the basic feature vector C = [450.031, 0.02, 0.03, 0.999, 0.000, 0.000, 0.0027].

[0087] S2, Based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data, construct a dynamic spatiotemporal graph network G;

[0088] The purpose of step S2 is to establish a graph structure representation for slope point cloud data that can simultaneously encode spatial topological relationships and temporal evolutionary correlations, enabling subsequent graph neural network inference to perceive the local geometric discontinuities and temporal evolutionary features of each point cloud node within a unified framework. Traditional deep learning methods (such as convolutional neural networks) rely on regular gridded data structures, while 3D point clouds are essentially unordered discrete point sets, unsuitable for directly applying standard convolution operations. By constructing a dynamic spatiotemporal graph network with point cloud nodes as graph nodes, local normal vector angle consistency as spatial edge weights, and temporal correspondence as temporal edge weights, irregular point cloud data can be converted into a graph structured representation, allowing subsequent graph convolution operations to extract the local structural features and historical evolutionary features of each node while preserving spatial topology and temporal correlations.

[0089] Furthermore, by using the consistency of the local normal vector angle as the spatial edge weight, when the local normal vector directions of adjacent point cloud nodes change abruptly (such as the ground surfaces on both sides of a crack, where the local normal vector directions deflect each other), the consistency of the angle between them (inner product value cosθ) decreases, and the spatial edge weight decreases accordingly. This allows the graph network to naturally perceive the discontinuity of local geometry, which is an important geometric criterion for distinguishing between micro-fractured regions and continuous surfaces.

[0090] In this embodiment, a dynamic spatiotemporal graph network G is constructed based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data, including:

[0091] S21, obtain any two point cloud nodes i and j at the same time node, and calculate the inner product value cosθ of the local normal vectors of the two points;

[0092] S22, based on the inner product value cosθ and the Euclidean distance between any two point cloud nodes, calculate the spatial edge weight w. s ;

[0093] The above analysis, combining steps S21 and S22, comprehensively considers the consistency of the angle between the local normal vectors of adjacent nodes (using the inner product value) when determining the spatial edge weight. Quantization) and Euclidean distance The two factors are calculated using the following formula: ,in Let i be the spatial edge weights of two point cloud nodes i and j. Let be the local normal vector of node i. Let be the local normal vector of node j. It is the inner product of the local normal vectors of the two points (within the range [-1, 1]). Normalize it to [0,1]. The Euclidean distance between the two points is in meters, and ε is a small positive number to prevent division by zero (typically ε = 0.01m), read from a preset database; this formula ensures that the local normal vector directions are highly consistent. →1) and the spatial edge weights between closely spaced adjacent nodes are higher, reflecting that they are in a geometrically continuous surface region, while the local normal vector direction changes abruptly ( →-1, such as between two sides of a crack) or between nodes that are far apart, have a lower spatial edge weight.

[0094] S23, connect point cloud nodes representing the same spatial physical location at adjacent time nodes to generate time edge weights w. t ;

[0095] Specifically, the temporal edge weights are generated by connecting point cloud nodes representing the same spatial physical location at different time points (the correspondence is determined by nearest neighbor matching). The normalized historical displacement rate amplitude is used as the temporal edge weight to encode the evolution intensity of the nodes in time.

[0096] For example, in the slope monitoring scenario described above, the local normal vector of node P (whose basic feature vector has been calculated in step S1) is: The local normal vector of its neighboring node Q is (The local slope normal vector shows a significant deflection, possibly located at the edge of a crack), Euclidean distance between the two nodes. m, then The spatial edge weights of nodes P and Q are represented by the following formula: The inner product value ,therefore In contrast, if the local normal vector of another neighboring node R of node P is If the distance is the same (0.35m), then The weights of the spatial edges between nodes P and R are represented by the following formula: Slightly higher This reflects the modulating effect of geometric discontinuity at the crack edge on the spatial edge weight.

[0097] S24, taking the point cloud nodes as node feature inputs, combined with spatial edge weights w s and time edge weight w t Establish a global spatiotemporal graph matrix to generate a dynamic spatiotemporal graph network G.

[0098] Specifically, the steps for generating the dynamic spatiotemporal graph network G are as follows:

[0099] Step A: Construct the multidimensional node feature matrix X

[0100] Extract a preset number of continuous time nodes ( Given the total number of consecutive time nodes, all point cloud nodes participating in the construction of the graph network (assuming the total number of point cloud nodes in a single period is N); arrange the basic feature vector corresponding to each point cloud node according to spatial index and time series, generating a dimension of The multidimensional node feature matrix X, where F is the length of the basic feature vector C; this multidimensional node feature matrix X constitutes the basic data input stream of the dynamic spatiotemporal graph network.

[0101] Step B: Construct the spatial adjacency matrix

[0102] For any independent time node t (where Using all point cloud nodes at that time point as a reference, initialize a dimension as follows: Spatial adjacency matrix .

[0103] Iterate through all node pairs at this time point. If there is a spatial edge weight generated by calculation between the i-th point cloud node and the j-th point cloud node... Then, the weight value is assigned to the corresponding element in the spatial adjacency submatrix, that is, let... If two nodes do not belong to the preset neighborhood in space, the corresponding element is assigned a value of 0. Thus, the generation... A set of independent spatial adjacency submatrices is used to represent the spatial topological connectivity of the slope surface at each time point in the graph network.

[0104] Step C: Construct the temporal adjacency submatrix

[0105] For two adjacent time points t and t+1, the system initializes a dimension as follows: Temporal adjacency matrix .

[0106] Due to the time-side weight w t Aimed at tracking the evolution of the same spatial physical location at different time points, the system extracts the connection relationship between the i-th point cloud node at the same spatial physical location at time point t and time point t+1, and applies the pre-generated time edge weights corresponding to node i. Assign values ​​to the diagonal elements of the temporal adjacency submatrix, i.e., set... Then, the off-diagonal elements of the matrix are all assigned the value 0.

[0107] Step D: Assemble the global spatiotemporal graph matrix

[0108] After generating the spatial and temporal adjacency submatrices, the spatial and temporal adjacency submatrices are concatenated in blocks to construct a matrix with dimension [missing information]. Global spatiotemporal graph matrix ;

[0109] The specific assembly rules are as follows: [The following text appears to be incomplete and requires further context:] ( (Total number of consecutive time nodes) spatial adjacency submatrices Arranged sequentially on the diagonal blocks of the global spatiotemporal graph matrix; temporal adjacency submatrices The matrix and its transpose are arranged on the second diagonal blocks on both sides of the diagonal block; the remaining irrelevant blocks are padded with 0.

[0110] Global Spatiotemporal Graph Matrix The block matrix form is as follows:

[0111] ;

[0112] This global spatiotemporal graph matrix fully records all boundary weight constraints of slope deformation in three-dimensional space and one-dimensional time.

[0113] Step E: Generate a dynamic spatiotemporal graph network G

[0114] The multidimensional node feature matrix X and the global spatiotemporal graph matrix Binding and encapsulation are performed to generate a dynamic spatiotemporal graph network. .

[0115] In computer memory, a dynamic spatiotemporal graph network G is strictly defined as an ordered set of data tuples. .

[0116] S3, based on the dynamic spatiotemporal graph network G, performs spatiotemporal feature fusion on the basic feature vector C to obtain the structural evolution index E of each point cloud node;

[0117] In this embodiment, the structural evolution index E refers to a numerical parameter that quantifies the degree of joint spatial anisotropy and temporal monotonicity of each point cloud node. The numerical value is constrained within a preset numerical range (such as [0,1]). This index is the core discriminant of the present invention for distinguishing micro-fracture features from vegetation noise.

[0118] The above S3 uses spatial message passing and time series analysis of graph neural networks to extract the structural evolution index E from the basic feature vector C, which can jointly represent "spatial anisotropy" and "temporal monotonicity". This index serves as the core driving force for the subsequent adaptive filtering mechanism. By fusing these two features through normalized multiplication, the obtained structural evolution index helps to distinguish between "micro-fracture evolution signals with spatiotemporal structural regularities" and "spatiotemporally randomly distributed vegetation noise", providing accurate node-level quantitative discrimination basis for subsequent adaptive filtering.

[0119] In this embodiment, based on the dynamic spatiotemporal graph network G, spatiotemporal feature fusion is performed on the basic feature vector C to obtain the structural evolution index of each point cloud node, including:

[0120] S31 uses the spatial convolutional layer in the dynamic spatiotemporal graph network and performs neighborhood node feature aggregation on the basic feature vector according to the spatial edge weights to extract the spatial anisotropic feature data of the generated point cloud nodes.

[0121] Specifically, spatial anisotropy refers to the degree of geometric asymmetry of the local surface around a point cloud node. At the edge of a crack or in a stress concentration area, there are obvious directional differences in the local normal vector direction of the surrounding neighboring points, and the degree of anisotropy is high; while on a flat slope or in an area with uniform vegetation distribution, the degree of anisotropy is low.

[0122] Specifically, the spatial convolutional layer aggregates the basic features of neighboring nodes through graph convolution operations, and the weighted aggregation formula is as follows: Where N(i) is the set of spatial neighbor nodes of node i. Let W be the spatial edge weights of two point cloud nodes i and j, and W be the learnable weight matrix. j Let σ be the basic feature vector of the neighboring node j, and σ be the activation function (such as ReLU). The output of this operation is... This refers to the spatial anisotropy characteristic data of node i. The larger the value, the more asymmetrical the geometry around node i is.

[0123] S32 extracts the temporal monotonic cumulative feature data of the generated point cloud nodes by performing temporal series gradient integration on the instantaneous displacement rate in the basic feature vector based on the temporal edge weights through the temporal series convolutional layer in the dynamic spatiotemporal graph network; wherein, the temporal series convolutional layer and the spatial convolutional layer are parallel branches.

[0124] Specifically, temporal monotonicity refers to whether the instantaneous displacement rate of a point cloud node shows a continuous increasing trend in the time dimension: the true micro-fracture precursor deformation usually shows a monotonic evolution that gradually accelerates over time, while vegetation noise shows an irregularity of random fluctuation in time series.

[0125] Specifically, the time-series convolutional layer performs gradient integration on the instantaneous displacement rate sequence of nodes at consecutive time points to quantify the cumulative growth trend of the instantaneous displacement rate: ,in: For the time-monotonically accumulated feature data of node i, The total number of consecutive time points. Let i be the time edge weight. Let Z be the Z-axis component of the instantaneous displacement velocity of node i at time node t. Let Z be the Z-axis component of the instantaneous displacement velocity of node i at time node t-1. This represents the rate difference between adjacent time steps; a positive value indicates that the rate is monotonically increasing.

[0126] S33 performs a normalized product operation on the spatial anisotropic feature data and the temporal monotonic cumulative feature data, and outputs a structural evolution index with numerical constraints within a preset numerical range; the structural evolution index is a joint data parameter characterizing the spatial anisotropy and temporal monotonic cumulativeity of point cloud nodes.

[0127] Specifically, the structural evolution index value of node i is obtained by performing a normalized product operation on the aforementioned two types of features. , ,in, For the spatial anisotropy feature data of node i, For the time-monotonically accumulated feature data of node i; Perform normalized product operation on two types of feature data. Let it be its global maximum value. By dividing by this global maximum value, the calculated result is... The structural evolution index value of node i is constrained within the interval [0,1]. For example, continuing the slope monitoring scenario above, node P, located at the edge of an expanding crack, has a spatial anisotropy characteristic normalized value of 0.85 (significant differences in the local normal vector directions on both sides of the crack) and a temporal monotonic cumulative characteristic normalized value of 0.78 (the Z-axis component of the instantaneous displacement rate in the past 5 periods was 0.8, 1.1, 1.5, 2.0, and 2.7 mm / day, respectively, showing a monotonic accelerating trend). , Let be the structural evolution index value of node P; while node Q, located in the shrub-covered area, has a spatial anisotropy normalized value of 0.20 (random vegetation distribution, weak directionality) and a temporal monotonic cumulative normalized value of 0.12 (random fluctuations in instantaneous displacement rate, no monotonicity). , The structural evolution index value of node Q is consistent with the physical characteristics of real vegetation noise, and the difference between the two is obvious, providing a clearer numerical basis for subsequent threshold judgment.

[0128] S4. Based on the structural evolution index, determine the data fidelity stiffness parameters in the preset local energy functional equation. ;

[0129] In this embodiment, the local energy functional equation consists of a data fidelity term weighted by data fidelity stiffness parameters and a spatial smoothing term weighted by preset smoothing weights.

[0130] Among them, data fidelity stiffness parameters It refers to the weight parameter in the local energy functional equation that measures the constraint strength of the original elevation data on the final filtering result. The larger the value, the closer the filtering result is to the original observation value; conversely, the smaller the value, the closer it is to the smooth result of the surrounding space.

[0131] Specifically, the general form of the local energy functional equation is: ,in The solution is the actual ground elevation data to be solved, and z is the original elevation data (Z coordinate of the point cloud node). To ensure data fidelity of stiffness parameters, To preset smoothing weights, For data fidelity items (measurement) (the degree of deviation from z) For spatial smoothness (measuring the magnitude of the spatial gradient of u, i.e., the surface roughness); when Much larger At the same time, optimize the driver. Approaching z, the filtering result highly preserves the original data; when much smaller At the same time, optimize the driver. It achieves deep noise reduction by approximating the smoothness of the surrounding space.

[0132] In this embodiment, the data fidelity stiffness parameters in the preset local energy functional equation are determined based on the structural evolution index, including:

[0133] If the structural evolution index is greater than the preset evolution threshold, the data fidelity stiffness parameter will be mapped to the preset maximum value and the smoothing exemption mechanism for micro-fracture data will be triggered; otherwise, the current structural evolution index will be mapped to the data fidelity stiffness parameter.

[0134] Since the calculation results of the structural evolution index alone cannot directly provide the filtered true surface elevation data, a mathematical framework is needed to combine this index with the surface elevation extraction process of point cloud data. The local energy functional equation is an optimization model under a variational framework. Its basic logic is to formulate the filtering problem as: seeking the optimal balance between ensuring the filtered result is faithful to the original observation data (data fidelity constraint) and ensuring the filtered result is spatially smooth (spatial smoothness constraint); this is achieved by mapping the structural evolution index output in step S3 to data fidelity stiffness parameters. It can assign differentiated data fidelity strength to point cloud nodes in different regions: for nodes with obvious structural evolution (high E value), it can improve... This ensures the filtered results closely match the original observations, thus preserving micro-fracture characteristics; for nodes with weak structural evolution (low E values), the filtering effect is reduced. The spatial smoothing term is allowed to play a dominant role, achieving effective removal of vegetation noise. This node-by-node adaptive stiffness control strategy enables the filter of this application to simultaneously achieve "preservation of micro-fracture features" and "removal of vegetation noise" in the same processing process, breaking through the fundamental limitations of traditional static filters.

[0135] In this embodiment, the current structural evolution index is mapped to a data-fidelity stiffness parameter. The specific process is as follows:

[0136] The system reads preset basic stiffness data and preset proportional constants from the preset database, which are usually set and stored in the preset database by the preset staff; inputs the structural evolution index into the system's preset tangent function model to calculate and output evolution increment data; multiplies the preset proportional constants with the evolution increment data to obtain the dynamic stiffness value; adds the dynamic stiffness value to the preset basic stiffness data to generate data-fidelity stiffness parameters.

[0137] It should be noted that the formula for calculating the data fidelity stiffness parameter is: ,in The preset basic stiffness data is read from the preset database, k0 is the preset proportional constant (read from the preset database), and E is the structural evolution index (value range [0,1]). The output of the tangent function model, i.e., the evolutionary increment data, is characterized by monotonically increasing from 0 to tan(π / 4)=1.0 within the interval E∈[0,1]. It exhibits a nonlinear amplification effect on E values ​​close to 1, making the stiffness increment in high evolutionary index regions more pronounced. For example, in the slope monitoring scenario described above, let... =1.0, k0=9.0, for node , The structural evolution index value for node P: , For node P, the data-fidelity stiffness parameter is provided; for node P... , The structural evolution index value of node Q; Provide the data-fidelity stiffness parameters for node Q;

[0138] The data fidelity stiffness of node P is approximately 5.3 times that of node Q, enabling differentiated processing of micro-fractured regions (high stiffness, original data is preserved) and vegetated areas (low stiffness, deep smoothing and noise reduction).

[0139] In this embodiment, mapping the data fidelity stiffness parameter to a preset maximum value includes:

[0140] S41, read the preset maximum value stored in the preset database. The magnitude of the preset maximum value is greater than the preset smoothing weight corresponding to the spatial smoothing term in the local energy functional equation. It is generally preset by the staff in advance.

[0141] S42, If the structural evolution index is greater than the preset evolution threshold, then perform the data assignment operation and allocate the preset maximum value to the data fidelity stiffness parameter.

[0142] S43, substitute the assigned data fidelity stiffness parameter into the data fidelity term of the local energy functional equation, and calculate the elevation residual data between the real surface elevation data to be solved and the original elevation data in each iteration of the preset solution algorithm;

[0143] S44. Numerical penalty calculation is applied to the elevation residual data according to the assigned data fidelity stiffness parameter. When the energy value of the local energy functional equation converges to the minimum value, the value of the output real surface elevation data is directly constrained to be equal to the original elevation data. Thus, the original elevation data that has not undergone smooth variation is retained and output as exempted micro-fracture characteristic data.

[0144] The above steps detail the mathematical implementation logic of the smoothing exemption mechanism. Its core is to apply a numerical penalty to the elevation residual data using a preset maximum value M, thereby achieving a hard constraint on the original elevation data within the energy minimization framework. Let the preset maximum value be... Preset smoothing weights =1.0, when the structural evolution index value of a certain node i is 1.0. When the preset evolution threshold is exceeded, the data fidelity stiffness parameter of node i is assigned a value. Substituting this maximum stiffness into the optimization equation, the optimal solution of the objective functional is: (when (when), the error is less than ,in, The optimal solution of the objective functional with respect to node i is the final output of the true surface elevation data. The original elevation data for node i; The neighborhood smoothing term represents the combined influence of neighboring nodes on the smoothing constraints of the current node i. It is typically obtained by weighted summation of the true elevation values ​​of the neighboring nodes. The neighborhood smoothing coefficient represents the weight coefficient related to the current node in the spatial smoothing term. It is typically the sum of the preset smoothing weight λ and the neighborhood weights, achieving a hard constraint effect where the numerical accuracy of the real surface elevation data is equal to that of the original elevation data. The preset evolution threshold determines the trigger boundary of the exemption mechanism. If the threshold is set too high, some real micro-fracture areas may be smoothed because their E-values ​​do not reach the threshold; if the threshold is set too low, vegetation noise areas may be misjudged as micro-fracture areas and thus exempted from smoothing. The preset evolution threshold can be calibrated by analyzing historical monitoring data with or without known micro-fracture characteristics. Typically, the 95th percentile of the E-value distribution in the vegetation area of ​​the E-value distribution histogram is used as a reference for the threshold to control the misjudgment rate within a reasonable range. In actual engineering, an initial value (e.g., 0.5) can be set first and then manually corrected based on the actual monitoring results. For example, continuing the above slope monitoring scenario, assuming the preset evolution threshold is 0.5, the structural evolution index value of node P... This triggers the exemption mechanism, affecting the data fidelity stiffness parameters of node P. The final output is the actual surface elevation data of node P. ≈450.031m (identical to the original elevation data), its millimeter-level bulge features were preserved without loss; the structural evolution index value of node Q Without triggering the exemption mechanism, the data-fidelity stiffness parameters of node Q are calculated using the tangent function model in step S4. After solving the energy minimization problem, the elevation contribution of the vegetation canopy in the original elevation is effectively smoothed out, and the true surface elevation data is output.

[0145] Furthermore, after mapping the data fidelity stiffness parameter to a preset maximum value, it also includes:

[0146] If the structural evolution index is less than the preset denoising threshold, and the preset denoising threshold is less than the preset evolution threshold, then the data fidelity stiffness parameter is mapped to the preset minimum value, so that the real surface elevation data is close to the macroscopic low-frequency elevation data output by the spatial smoothing term, triggering a deep denoising mechanism for vegetation cover data.

[0147] If the structural evolution index is between the preset denoising threshold and the preset evolution threshold, then linear interpolation is performed based on the current value of the structural evolution index, and the output is a data-fidelity stiffness parameter with a linear transition.

[0148] The above steps further extend the adaptive mapping mechanism of data fidelity stiffness parameters, expanding the single threshold judgment into a three-interval hierarchical processing strategy, forming a complete three-stage control logic of "preservation-transition-denoising". When the structural evolution index value of node i is greater than the preset evolution threshold (first interval), the smoothing exemption mechanism is triggered. (Preset maximum value), where For the data fidelity stiffness parameter of node i; when When the preset denoising threshold is reached (second interval), the deep denoising mechanism is triggered. ( To minimize the noise, such as m=0.001, the actual surface elevation data is made closer to the macroscopic low-frequency elevation data (i.e., spatial smoothing result), and high-frequency noise from vegetation is deeply removed; when the preset noise reduction threshold is ≤ When the value is less than or equal to the preset evolution threshold (third interval, transition zone), the stiffness parameters for linear transition are output through linear interpolation: ,in, For the data fidelity stiffness parameters of node i, Let i be the structural evolution index value of node i. To preset the noise reduction threshold, To preset the evolution threshold, and The values ​​of the data fidelity stiffness parameters corresponding to the two thresholds are respectively; the method for determining the preset denoising threshold is: the mean of the vegetation area in the E-value distribution histogram plus one standard deviation, usually between 0.05 and 0.15, to ensure that the E-value of the vegetation-covered area is generally lower than this threshold; this three-stage control strategy avoids the "abrupt boundary" effect caused by a single threshold judgment, making the spatial distribution of stiffness parameters smoother and more continuous, thereby improving the physical rationality of the filtering results.

[0149] In one specific embodiment, for example, let the preset denoising threshold be 0.10 and the preset evolution threshold be 0.50. =0.001, =6.0, then for =0.30 (in the transition zone): This reflects the smooth transition characteristic of moderate structural evolution corresponding to moderate data fidelity; if =0.05 (most likely a densely vegetated area): =0.001, the actual surface elevation data output will be closer to the spatially smooth macroscopic low-frequency elevation, effectively removing the vegetation canopy elevation deviation.

[0150] S5 updates the local energy functional equation based on the data-fidelity stiffness parameter, solves the minimum value of the updated local energy functional equation through a preset solution algorithm, and outputs real surface elevation data covering the smoothed area and the exempted retention area.

[0151] Step S5 is an important supplement and enhancement to the adaptive stiffness control mechanism. The data fidelity stiffness parameter is smoothly mapped through the tangent function model. Its maximum value is limited by the numerical range of the tan function and the preset proportional constant. For nodes whose structural evolution indicators have exceeded the preset evolution threshold, the data fidelity provided by this smoothing mapping may still be insufficient to completely suppress the modification of the original elevation data by the spatial smoothing term, thus causing extremely subtle micro-fracture features to be partially smoothed. Therefore, the data fidelity stiffness parameter of the node exceeding the threshold is directly mapped to the preset maximum value (whose numerical magnitude is far greater than the smoothing weight of the spatial smoothing term) through the smoothing exemption mechanism. Mathematically, the influence of the spatial smoothing term is completely suppressed, so that the final output real surface elevation data of the node is forced to be equal to the original elevation data in the energy minimization iteration, thereby achieving lossless preservation of micro-fracture feature data.

[0152] In this embodiment, the minimum value of the updated local energy functional equation is solved using a preset solution algorithm, including:

[0153] S51, Substitute the original elevation data and data fidelity stiffness parameters into the data fidelity term to establish the first polynomial describing the data fidelity error;

[0154] S52, Substitute the real surface elevation data to be solved and the spatial neighborhood elevation data of the point cloud nodes into the spatial smoothing term to establish a second polynomial describing the surface smoothing error.

[0155] S53, combine the first polynomial and the second polynomial to transform them into the objective functional matrix;

[0156] S54. Based on the preset sparse matrix optimization algorithm, the target functional matrix is ​​differentiated and iteratively degraded. When the iteration residual is less than the preset convergence threshold, the real surface elevation data that meets the energy minimization condition is calculated.

[0157] The above steps provide a detailed solution process for finding the local energy functional equation minimum. The core of this process is transforming the variational optimization problem into a numerical solution of a sparse linear equation system. In a specific embodiment, let the total number of point cloud nodes be N, which is the total number of point cloud nodes in a single period. The actual surface elevation data of each node are then arranged into a column vector. , These represent the actual surface elevation data from the 1st to the Nth data points to be solved, where T is the transpose operator, and the original elevation data are arranged as follows: , These represent the original elevation data from the 1st to the Nth time, with the data fidelity stiffness parameters arranged in a diagonal matrix. Let represent the data fidelity stiffness parameters from the 1st to the Nth, respectively. Then, the matrix form of the objective functional equation is: ,in The solution is the actual ground elevation data to be solved, and z is the original elevation data (Z coordinate of the point cloud node). To predetermine smoothing weights, L is the graph Laplacian matrix, representing the difference in elevation between each node and its spatial neighbors; differentiate the above objective functional equation and let... This yields the corresponding system of linear equations: The coefficient matrix in the above system of equations For sparse symmetric positive definite matrices (where non-zero elements appear only at the corresponding positions of adjacent nodes), efficient solutions can be obtained using sparse matrix optimization algorithms such as the conjugate gradient (CG) method or the preconditional conjugate gradient method; the convergence threshold is usually set to the relative residual. ( For the first The residual vector of the next iteration (r0 is the initial residual vector of the first iteration). When this condition is met, the iteration stops, and the true surface elevation data that satisfies the energy minimization condition is output. .

[0158] For example, in a local region containing 1000 nodes, the coefficient matrix (A+λL) is a sparse matrix of 1000×1000, where the number of non-zero elements is approximately 5 to 10 times the number of nodes (depending on the average number of neighborhood points). The conjugate gradient method typically converges to a result that meets the accuracy requirements within 50 to 200 iterations; for node P ( , m, the neighborhood smoothed elevation mean is approximately 449.70m, λ=1.0), its true surface elevation data The bulge feature at 31 mm in elevation direction was highly preserved, with only a minor correction of about 21 mm due to neighborhood constraints.

[0159] S6, based on multiple sets of real surface elevation data, reconstructs high-fidelity digital elevation models corresponding to each time node, and performs data spatial difference calculation on the high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data.

[0160] In this embodiment, the high-fidelity digital elevation model (DEM) refers to pure surface elevation raster data that has been processed by adaptive filtering to retain micro-fracture features while removing interference from non-ground elements such as vegetation. It is the basic data for subsequent deformation calculations.

[0161] Step S6 converts the real surface elevation data for each period, after adaptive filtering, into pure three-dimensional displacement field data that can be quantitatively analyzed, providing a quantitative basis for the final judgment of slope instability. Since steps S4 to S5 have completed the core processing of removing vegetation noise while preserving micro-fracture characteristics, the high-fidelity digital elevation model reconstructed based on the real surface elevation data more accurately represents the pure surface morphology of each period, in which the millimeter-level bulges and cracks in the micro-fracture areas are realistically reflected. By performing spatial difference calculations on the high-fidelity digital elevation models of adjacent periods, the elevation change of each spatial location between the two periods can be quantified more accurately, thereby constructing pure three-dimensional displacement field data containing spatial coordinates and vertical deformation information.

[0162] Specifically, spatial difference calculations are performed on high-fidelity digital elevation models at adjacent time points to generate clean three-dimensional displacement field data, including:

[0163] S61, extract the first digital elevation model data corresponding to the first time node from the pre-stored high-fidelity digital elevation model, and extract the second digital elevation model data corresponding to the second time node adjacent to the first time node;

[0164] S62, according to the preset spatial grid resolution, perform spatial grid projection and alignment processing on the first digital elevation model data and the second digital elevation model data to generate aligned grid data with the same plane coordinate sequence;

[0165] S63, loop through each plane coordinate point in the aligned grid data, and extract the first elevation data of the plane coordinate point in the first digital elevation model data and the second elevation data in the second digital elevation model data respectively;

[0166] S64, perform a numerical subtraction operation between the second elevation data corresponding to the same plane coordinate point and the first elevation data to calculate and obtain the elevation change data corresponding to each plane coordinate point;

[0167] It should be noted that the formula for elevation change data is: ,in The first time node (time node) The first elevation data at the plane coordinates (x, y). The second time node (time node) , The second elevation data at the same plane coordinates, This represents the elevation change data at coordinates (x, y). Positive values ​​indicate an increase in surface elevation (such as bulging), while negative values ​​indicate a decrease in elevation (such as subsidence or erosion).

[0168] S65 extracts the calculated elevation change data and performs a three-dimensional vector splicing operation on the elevation change data and the corresponding plane coordinate points to generate pure three-dimensional displacement field data containing spatial coordinate information and vertical deformation values.

[0169] It should be noted that the formula for constructing the pure three-dimensional displacement field data D(x,y) is: , where T is the transpose operator; spatial grid alignment uses bilinear interpolation to resample the two high-fidelity digital elevation models onto a unified planar coordinate grid (the preset spatial grid resolution is usually set to 0.1m×0.1m or 0.05m×0.05m, and the choice of resolution should take into account the UAV point cloud density and computing resources), ensuring that the subtraction operation is performed at strictly corresponding spatial locations.

[0170] Continuing with the aforementioned slope monitoring scenario, if the area of ​​the target slope monitoring region is... With a spatial grid resolution of 0.1m, a total of 1000 × 500 = 500,000 planar coordinate points are generated. At the coordinate point (100.0m, 200.0m), the elevation h1 of the high-fidelity digital elevation model in phase 4 is 450.010m, and the elevation h2 of the high-fidelity digital elevation model in phase 5 is 450.052m. Therefore, the elevation change data Δh = +0.042m (+42mm) indicates that a surface uplift of 42mm occurred at this location between phases 4 and 5. Considering this location... With a preset evolution threshold, the uplift is judged as a genuine micro-fracture bulge feature rather than vegetation noise. Finally, pure three-dimensional displacement field data D=[100.0,200.0,0.042]m is generated by splicing, providing more accurate deformation input for subsequent instability judgment.

[0171] S7. In the pure three-dimensional displacement field data, if the target deformation rate in the region where the structural evolution index is greater than the preset evolution threshold exceeds the preset instability critical value, then the slope instability early warning data will be output to the terminal.

[0172] Step S7 above only triggers an early warning in the intersection area that meets two conditions: First, the structural evolution index E of the area exceeds the preset evolution threshold, proving that the area is undergoing genuine structural evolution (rather than vegetation noise interference); Second, the target deformation rate of the area exceeds the preset instability critical value, proving that the deformation rate has reached a dangerous level indicating that the slope is about to become unstable; whereby the formula for calculating the target deformation rate is: ,in This represents the elevation change data at coordinates (x, y), where Δt is the time interval between two adjacent time nodes. The preset instability critical value should be determined comprehensively based on engineering geological parameters such as slope type (soil slope, rock slope), slope angle, and hydrological conditions. It can be set with reference to relevant technical specifications or critical deformation rate data from similar historical cases. Generally, the critical deformation rate of soil slopes is lower than that of rock slopes. The specific value can be calibrated in conjunction with the actual situation in engineering applications. When the target deformation rate in the area where the structural evolution index exceeds the threshold continues to increase but has not yet exceeded the preset instability critical value, warning level data (such as yellow, orange, and red graded warnings) can be further output so that engineering management personnel can take graded response measures.

[0173] For example, continuing the slope monitoring scenario described above, assuming a preset instability threshold of 5 mm / day, and a time interval of 7 days between two adjacent time points, at coordinates (100.0m, 200.0m)... If the elevation change data Δh = +42mm in the 4th to 5th periods (preset evolution threshold), then the target deformation rate at coordinates (100.0m, 200.0m) = 42 / 7 = 6mm / day > 5mm / day (preset instability threshold), triggering a slope instability early warning output. The early warning data (including the coordinates of the early warning location, the value of the target deformation rate, the value of the structural evolution index, etc.) is sent to the monitoring terminal, prompting the engineering management personnel that there is a risk of slope instability at this location and that emergency measures must be taken immediately.

[0174] Example 2

[0175] like Figure 3 As shown, the difference between this embodiment and Embodiment 1 is that this embodiment provides a slope instability early warning system based on UAV remote sensing data. This system corresponds one-to-one with the slope instability early warning method based on UAV remote sensing data in Embodiment 1. The system includes:

[0176] The data acquisition module is used to acquire three-dimensional point cloud data of the target slope after spatial registration at a preset number of continuous time nodes. The three-dimensional point cloud data includes the local normal vector of each point cloud node and the basic feature vector of the original elevation data.

[0177] The feature processing module is used to construct a dynamic spatiotemporal graph network based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data; based on the dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vectors to obtain the structural evolution index of each point cloud node.

[0178] The filtering control module is used to determine the data fidelity stiffness parameters in the preset local energy functional equation based on the structural evolution index; update the local energy functional equation based on the data fidelity stiffness parameters; solve the minimum value of the updated local energy functional equation through a preset solution algorithm; and output the real surface elevation data covering the smoothed area and the exempted retention area.

[0179] The early warning output module is used to reconstruct high-fidelity digital elevation models corresponding to each time node based on multiple sets of real surface elevation data, and to perform data spatial difference calculation on the high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data. In the pure three-dimensional displacement field data, if the target deformation rate in the area where the structural evolution index is greater than the preset evolution threshold exceeds the preset instability critical value, the slope instability early warning data is output to the terminal.

[0180] The execution process of each unit can be carried out according to the steps of the slope instability early warning method based on UAV remote sensing data in Example 1, and will not be described in detail in this example.

[0181] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0182] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0183] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0184] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0185] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A slope instability early warning method based on UAV remote sensing data, characterized in that, The method includes: Acquire spatially registered three-dimensional point cloud data of the target slope at a preset number of continuous time nodes. The three-dimensional point cloud data includes the local normal vector of each point cloud node and the basic feature vector of the original elevation data. Based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data, a dynamic spatiotemporal graph network is constructed. Based on the dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vector to obtain the structural evolution index of each point cloud node. Based on the structural evolution index, determine the data fidelity stiffness parameters in the preset local energy functional equation; The local energy functional equation is updated based on the data-fidelity stiffness parameter. The minimum value of the updated local energy functional equation is solved by a preset solution algorithm, and the output is real surface elevation data covering the smoothed area and the exempted retention area. Based on multiple sets of real surface elevation data, high-fidelity digital elevation models corresponding to each time node are reconstructed, and data spatial difference calculations are performed on the high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data. In pure three-dimensional displacement field data, if the target deformation rate in a region where the structural evolution index is greater than the preset evolution threshold exceeds the preset instability critical value, slope instability early warning data will be output to the terminal.

2. The slope instability early warning method based on UAV remote sensing data according to claim 1, characterized in that, The steps for obtaining the local normal vectors of each point cloud node and the basic feature vectors of the original elevation data are as follows: At the current time point, extract the three-dimensional spatial coordinates of the point cloud nodes and generate neighboring point set data based on the preset neighborhood radius; A local covariance matrix is ​​constructed based on the three-dimensional spatial coordinates of the adjacent point set data. The local covariance matrix is ​​then decomposed into eigenvalues, and the eigenvectors are extracted as local normal vectors. Extract the difference in three-dimensional spatial coordinates of point cloud nodes between the current time node and the previous time node, and divide the difference in three-dimensional spatial coordinates by the corresponding time interval data to obtain the instantaneous displacement rate; The basic feature vector is generated by splicing together the original elevation data corresponding to the three-dimensional spatial coordinates, the local normal vector data, and the instantaneous displacement rate.

3. The slope instability early warning method based on UAV remote sensing data according to claim 1, characterized in that, Based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data, a dynamic spatiotemporal graph network is constructed, including: Get any two point cloud nodes at the same time point and calculate the inner product of the two point cloud nodes. Based on the inner product value and the Euclidean distance data between any two point cloud nodes, calculate and generate spatial edge weights; Connect point cloud nodes representing the same spatial physical location at adjacent time points to generate time edge weights; The point cloud nodes are used as node feature inputs, and a global spatiotemporal graph matrix is ​​established by combining the spatial edge weights and the temporal edge weights to generate a dynamic spatiotemporal graph network.

4. The slope instability early warning method based on UAV remote sensing data according to claim 3, characterized in that, The specific steps for generating a dynamic spatiotemporal graph network are as follows: Extract a preset number of continuous time nodes Down, The total number of consecutive time nodes is N, representing all point cloud nodes participating in the construction of the graph network. The total number of point cloud nodes in a single period is N. The basic feature vector C corresponding to each point cloud node is arranged according to spatial index and time series, generating a dimension of N. The multidimensional node feature matrix X; where F is the length of the basic feature vector C; Traverse all node pairs at any given time point and construct a spatial adjacency submatrix based on the spatial edge weights of the node pairs; the spatial adjacency submatrix is ​​used to represent the spatial topological connectivity of the slope surface at each time point in the graph network; Based on the connection relationship between any point cloud node at the same spatial physical location at time node t and time node t+1, construct a time adjacency submatrix; The spatial adjacency submatrix and the temporal adjacency submatrix are segmented and concatenated to construct a global spatiotemporal graph matrix. ; The multidimensional node feature matrix X and the global spatiotemporal graph matrix are combined. Binding and encapsulation are performed to generate a dynamic spatiotemporal graph network. .

5. The slope instability early warning method based on UAV remote sensing data according to claim 3, characterized in that, Based on a dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vectors to obtain the structural evolution indices of each point cloud node, including: The spatial anisotropic feature data of the generated point cloud nodes are extracted by using the spatial convolutional layer in the dynamic spatiotemporal graph network and by aggregating the neighborhood node features of the basic feature vector according to the spatial edge weights. By using the time-series convolutional layer in the dynamic spatiotemporal graph network and performing time-series gradient integration on the instantaneous displacement rate in the basic feature vector based on the time edge weights, the time monotonic cumulative feature data of the generated point cloud nodes is extracted. The spatial anisotropic feature data and the temporal monotonic cumulative feature data are subjected to normalized multiplication to output a structural evolution index with numerical constraints within a preset numerical range.

6. The slope instability early warning method based on UAV remote sensing data according to claim 1, characterized in that, Based on structural evolution indices, determine the data fidelity stiffness parameters in the preset local energy functional equations, including: If the structural evolution index is greater than the preset evolution threshold, the data fidelity stiffness parameter will be mapped to the preset maximum value and the smoothing exemption mechanism for micro-fracture data will be triggered; otherwise, the current structural evolution index will be mapped to the data fidelity stiffness parameter. The step of mapping the data fidelity stiffness parameter to a preset maximum value includes: Read the preset maximum value stored in the preset database. The magnitude of the preset maximum value is greater than the preset smoothing weight corresponding to the spatial smoothing term in the local energy functional equation. If the structural evolution index is greater than the preset evolution threshold, a data assignment operation is performed to allocate the preset maximum value to the data fidelity stiffness parameter. The assigned data fidelity stiffness parameter is substituted into the data fidelity term of the local energy functional equation, and in each iteration of the preset solution algorithm, the elevation residual data between the real surface elevation data to be solved and the original elevation data is calculated. Numerical penalty calculations are applied to the elevation residual data based on the assigned data fidelity stiffness parameters, and when the energy value of the local energy functional equation converges to a minimum, the output value of the true surface elevation data is directly constrained to be equal to the original elevation data.

7. The slope instability early warning method based on UAV remote sensing data according to claim 6, characterized in that, After mapping the data fidelity stiffness parameter to a preset maximum value, the following is also included: If the structural evolution index is less than the preset denoising threshold, and the preset denoising threshold is less than the preset evolution threshold, then the data fidelity stiffness parameter is mapped to the preset minimum value, so that the real surface elevation data approaches the macroscopic low-frequency elevation data output by the spatial smoothing term, triggering a deep denoising mechanism for vegetation cover area data. If the structural evolution index is between a preset denoising threshold and a preset evolution threshold, then linear interpolation is performed based on the current value of the structural evolution index to output a data fidelity stiffness parameter with a linear transition.

8. The slope instability early warning method based on UAV remote sensing data according to claim 1, characterized in that, The updated local energy functional equations are solved using a pre-defined solution algorithm, including: Substitute the original elevation data and data fidelity stiffness parameters into the data fidelity term to establish the first polynomial describing the data fidelity error. Substitute the actual surface elevation data to be solved and the spatial neighborhood elevation data of the point cloud nodes into the spatial smoothing term to establish a second polynomial describing the surface smoothing error. The first polynomial and the second polynomial are combined to transform the target functional matrix. The target functional matrix is ​​differentiated and iteratively degraded based on a preset sparse matrix optimization algorithm. When the iteration residual is less than a preset convergence threshold, the real surface elevation data that satisfies the energy minimization condition is calculated.

9. The slope instability early warning method based on UAV remote sensing data according to claim 1, characterized in that, Spatial difference calculations are performed on high-fidelity digital elevation models at adjacent time points to generate clean three-dimensional displacement field data, including: Extract the first digital elevation model data corresponding to the first time node from the pre-stored high-fidelity digital elevation model, and extract the second digital elevation model data corresponding to the second time node adjacent to the first time node; Based on the preset spatial grid resolution, the first digital elevation model data and the second digital elevation model data are subjected to spatial grid projection and alignment processing to generate aligned grid data with the same plane coordinate sequence. Iterate through each plane coordinate point in the aligned grid data, and extract the first elevation data of the plane coordinate point in the first digital elevation model data and the second elevation data in the second digital elevation model data respectively. The elevation change data for each plane coordinate point is calculated by subtracting the second elevation data from the first elevation data. The calculated elevation change data is extracted, and the elevation change data is combined with the corresponding plane coordinate points to perform a three-dimensional vector stitching operation to generate pure three-dimensional displacement field data containing spatial coordinate information and vertical deformation values.

10. A slope instability early warning system based on UAV remote sensing data, characterized in that, The system includes: The data acquisition module is used to acquire three-dimensional point cloud data of the target slope after spatial registration at a preset number of continuous time nodes. The three-dimensional point cloud data includes the local normal vector of each point cloud node and the basic feature vector of the original elevation data. The feature processing module is used to construct a dynamic spatiotemporal graph network based on the consistency of the angle between the local normal vectors of adjacent point cloud nodes in the 3D point cloud data; based on the dynamic spatiotemporal graph network, spatiotemporal feature fusion is performed on the basic feature vectors to obtain the structural evolution index of each point cloud node. The filtering control module is used to determine the data fidelity stiffness parameters in the preset local energy functional equation based on the structural evolution index; update the local energy functional equation based on the data fidelity stiffness parameters; solve the minimum value of the updated local energy functional equation through a preset solution algorithm; and output the real surface elevation data covering the smoothed area and the exempted retention area. The early warning output module is used to reconstruct high-fidelity digital elevation models corresponding to each time node based on multiple sets of real surface elevation data, and to perform data spatial difference calculation on the high-fidelity digital elevation models of adjacent time nodes to generate pure three-dimensional displacement field data. In the pure three-dimensional displacement field data, if the target deformation rate in the area where the structural evolution index is greater than the preset evolution threshold exceeds the preset instability critical value, the slope instability early warning data is output to the terminal.

Images