Slope displacement monitoring data processing system based on unmanned aerial vehicle laser radar

Abstract

The present application provides a slope displacement monitoring data processing system based on unmanned aerial vehicle laser radar, relates to the technical field of data processing, and comprises: a spatial interpolation processing device for processing terrain feature data sets and constructing a digital terrain surface model; a device for calculating a slope surface displacement vector field, identifying a potential sliding surface and a deformation anomaly area, and generating a displacement field solution result based on the digital terrain surface model through spatial similarity analysis and iterative optimization by fusing a gradient descent algorithm; and a device for inputting the displacement field solution result into a risk assessment model and performing quantitative evaluation of slope stability through a multi-source data fusion analysis platform. The present application improves the practicality and operability of monitoring results.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a slope displacement monitoring data processing system based on UAV lidar. Background Technology

[0002] For slopes on specific sections of mountain highways, slope stability monitoring is a core technical challenge to ensure project safety.

[0003] Existing data processing procedures lack a standardized system, and the monitoring accuracy and efficiency in complex environments may fail to meet engineering requirements. For example, current monitoring technologies rely on inclinometers or piezometers, which have the following limitations: the slope is long and high, and both inclinometers and piezometers are deployed at single points (usually monitoring sections are set up at certain intervals, with a small number of monitoring points at each section), only reflecting deep deformation and seepage pressure at local points. For instance, within a certain thickness of the lower silty clay, if a section of soil experiences a sudden increase in seepage pressure due to uneven distribution of boulders, a single-point piezometer may struggle to detect this anomaly, easily missing landslide warnings. Furthermore, the lower silty clay is loose and contains boulders, requiring drilling to deep layers for inclinometers. Drilling can damage the original soil structure, reducing local shear strength and potentially inducing small-scale landslides, contradicting the original purpose of slope monitoring and protection. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a slope displacement monitoring data processing system based on UAV lidar, which improves the practicality and operability of the monitoring results.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0006] The first aspect is a slope displacement monitoring data processing system based on UAV lidar, including:

[0007] The slope displacement monitoring data processing system based on UAV lidar includes:

[0008] The slope displacement monitoring data processing system based on UAV lidar includes:

[0009] The acquisition module is used to acquire the raw lidar observation dataset;

[0010] The processing module is used to preprocess the original lidar observation dataset to obtain an optimized lidar observation dataset.

[0011] The computation module extracts stable terrain feature points from the observation data. Based on curvature and normal vector change rate indices, it selects surface protrusions, rock strata boundaries, and abrupt topographic changes as high-reliability feature points, generating an initial feature point set. Feature descriptors are constructed from this initial set, generating feature description vectors based on local surface geometric properties, resulting in a standardized feature description set. The standardized feature description set undergoes similarity matching, and mismatched point pairs are eliminated through a bidirectional consistency check, generating a reliable feature point correspondence set. Based on this reliable feature point correspondence set, the final spatial transformation parameters are calculated, generating a coordinate transformation parameter set. The coordinate transformation parameter set is then used to rigidly transform all lidar observation data, uniformly converting them to a reference coordinate system, generating a spatially unified dataset. Finally, the spatially unified dataset is fused, eliminating redundant points in overlapping areas, generating a complete three-dimensional slope observation dataset.

[0012] The generation module is used to establish an directional reference frame based on a slope three-dimensional observation dataset, with two preset stable reference benchmarks, and to generate adaptive correction factors according to the topological relationships and structural characteristics of each region within the frame.

[0013] The classification module optimizes the land cover classification parameters based on an adaptive correction factor to achieve accurate separation of vegetation cover, artificial structures and natural land surface, and generate an optimized terrain feature dataset.

[0014] The building module is used to perform spatial interpolation on the terrain feature dataset to build a digital terrain surface model;

[0015] The analysis module is used to calculate the displacement vector field of the slope surface based on the digital terrain surface model, through spatial similarity analysis and gradient descent algorithm for iterative optimization, identify potential sliding surfaces and abnormal deformation areas, and generate displacement field solution results.

[0016] The output module is used to input the displacement field calculation results into the risk assessment model and perform quantitative evaluation of slope stability through a multi-source data fusion analysis platform.

[0017] Furthermore, the original lidar observation dataset is preprocessed to obtain an optimized lidar observation dataset, including:

[0018] The original lidar observation dataset is subjected to noise filtering, and discrete noise points are removed based on statistical outlier analysis to generate a preliminary clean dataset. The preliminary clean dataset is then subjected to intensity correction, and intensity distortion caused by measurement geometry effects is eliminated through distance and incident angle normalization to generate an intensity-normalized dataset.

[0019] Coordinate correction is performed on the intensity normalized dataset. Based on sensor pose parameters and inertial measurement unit data, geometric deviations caused by changes in the attitude of the UAV platform are compensated, and a geometric correction dataset is generated.

[0020] The geometric correction dataset is sampled and optimized, the point density is adaptively adjusted based on the slope topographic features, redundant data points are removed, and an optimized lidar observation dataset is generated.

[0021] Furthermore, based on the slope three-dimensional observation dataset, two stable reference benchmarks are preset, an directional reference frame is established, and adaptive correction factors are generated according to the topological relationships and structural characteristics of each region within the frame, including:

[0022] From the slope three-dimensional observation dataset, select the first spatial reference point located in the stable bedrock area at the top of the slope and the second spatial reference point located in the stable bedrock area at the toe of the slope to generate the initial reference axis.

[0023] Using the initial reference axis as the reference direction, a regional coordinate system is constructed with the first spatial reference point as the origin and the direction of the reference line as the main axis, thus generating an orientation reference frame;

[0024] Based on the directional reference frame, and according to the geological structure characteristics of the slope and the trend of topographic change, the monitoring area is divided into several independent analysis units, and the unit division results are generated.

[0025] Extract the terrain structure attributes of each analysis unit to generate a set of unit terrain feature descriptions; based on the relative spatial position relationship between each analysis unit and two spatial reference points, calculate the azimuth weight and distance weight of each unit relative to the reference line to generate a set of unit weight distributions.

[0026] By integrating the unit terrain feature description set and the unit weight distribution set, the comprehensive terrain stability index of each analysis unit is calculated, and a unit stability evaluation set is generated.

[0027] Based on the unit stability evaluation set, the stability index is converted into adjustment coefficients for land cover classification parameters, and an adaptive correction factor is generated.

[0028] Furthermore, the land cover classification parameters are optimized based on an adaptive correction factor to achieve accurate separation of vegetation cover, artificial structures, and natural surfaces, generating an optimized terrain feature dataset, including:

[0029] Based on the adaptive correction factor, the parameter adjustment coefficients of each analysis unit are analyzed to generate a set of unit parameter adjustments. Based on the set of unit parameter adjustments, the classification threshold parameters in the land cover classification algorithm are adaptively adjusted according to the terrain stability characteristics of different analysis units to generate an optimized set of classification parameters.

[0030] Preliminary land cover classification is performed on the slope 3D observation dataset based on the optimized classification parameter set to generate initial classification results; the initial classification results are then processed to generate optimized classification results.

[0031] Based on the optimized classification results, the spatial distribution information of vegetation cover area and artificial structure area is extracted to generate ground feature mask dataset; the ground feature mask dataset is used to process the slope three-dimensional observation dataset, and the data points corresponding to vegetation and artificial structures are removed to generate a three-dimensional point set of bare surface.

[0032] The three-dimensional point set of exposed ground surface is processed to generate an optimized terrain feature dataset.

[0033] Furthermore, the terrain feature dataset is spatially interpolated to construct a digital terrain surface model, including:

[0034] Topographic feature analysis is performed on the optimized topographic feature dataset. Topographic feature points are identified based on local curvature changes and elevation variation coefficients. Ridge lines, valley lines, and slope abrupt change points are selected as key topographic points to generate a candidate topographic feature point set. Density optimization sampling is performed on the candidate topographic feature point set. The sampling density is adaptively adjusted based on the topographic complexity, with sparse sampling in flat areas and dense sampling in complex areas to generate a structured topographic point set.

[0035] A spatial interpolation function is constructed based on a structured terrain point set to obtain the final interpolation parameters and generate an optimized interpolation function. The optimized interpolation function is then used to calculate the numerical values ​​of regularly distributed elevation points to generate a preliminary terrain surface model.

[0036] The preliminary terrain surface model is smoothed to eliminate local interpolation anomalies, generating an optimized terrain surface model. The optimized terrain surface model is then subjected to integrity verification, and data gaps are filled based on the principle of terrain continuity to generate the final digital terrain surface model.

[0037] Furthermore, based on the digital terrain surface model, spatial similarity analysis and gradient descent algorithm are used for iterative optimization to calculate the slope surface displacement vector field, identify potential sliding surfaces and abnormal deformation areas, and generate displacement field solution results, including:

[0038] The digital terrain surface model is divided into several analysis units. An elevation residual objective function is established in each analysis unit to generate a unit difference quantification dataset.

[0039] Based on the unit difference quantization dataset, the gradient descent algorithm is used for iterative optimization and solution. By calculating the gradient of the objective function and gradually adjusting the spatial position parameters of the analysis unit along the negative gradient direction, a set of unit displacement parameters is generated.

[0040] Based on the unit displacement parameter set, the displacement vectors of all analysis units are integrated to construct a global displacement vector field. After statistical consistency testing, an optimized displacement vector field is generated, and the spatial distribution characteristics of displacement are analyzed to generate potential sliding boundary identification results.

[0041] The results of potential sliding boundary identification are analyzed to determine the spatial morphology and deformation evolution characteristics of the sliding surface, and displacement field calculation results are generated.

[0042] Furthermore, the displacement field calculation results are input into the risk assessment model, and a quantitative evaluation of slope stability is conducted through a multi-source data fusion analysis platform, including:

[0043] Based on the displacement field calculation results, the displacement vector field data and potential sliding surface information contained therein are extracted to generate an initial risk assessment input set; based on the initial risk assessment input set, the geological survey data and geotechnical parameters of the slope area are integrated to generate a geomechanical parameter set;

[0044] By integrating geomechanical parameter sets with real-time hydrological and meteorological data, a multi-source fusion database is generated; based on the multi-source fusion database, the safety factor of each potential sliding surface is calculated, and preliminary stability assessment results are generated.

[0045] Based on the preliminary stability assessment results, stress-strain analysis results are generated; by integrating the stress-strain analysis results with the preliminary stability assessment results, a comprehensive stability index is calculated using a weighted fusion algorithm to generate a slope stability level evaluation set.

[0046] Based on the slope stability level evaluation set, slope risk areas are divided and early warning levels are determined, generating slope risk zoning and early warning schemes; based on the slope risk zoning and early warning schemes, a quantitative evaluation report on slope stability is generated.

[0047] Furthermore, based on the unit difference quantization dataset, an iterative optimization solution is performed using the gradient descent algorithm. By calculating the gradient of the objective function and gradually adjusting the spatial position parameters of the analysis unit along the negative gradient direction, a set of unit displacement parameters is generated, including:

[0048] Initialize the displacement parameters of each analysis unit to generate an initial displacement parameter set, and calculate the elevation residual objective function value of each analysis unit to generate an initial objective function value set;

[0049] The initial displacement parameter set is iteratively optimized. In each iteration, the partial derivative of the objective function of each analysis unit with respect to the displacement parameters is calculated, a gradient vector set is generated, and the update direction and step size of the displacement parameters of each analysis unit are determined. The displacement parameters are adjusted along the negative gradient direction to generate the updated displacement parameter set.

[0050] The updated displacement parameter set is recalculated to generate a new objective function value set; the new objective function value set is compared with the objective function value set of the previous iteration to determine whether the iteration meets the preset convergence condition, and finally the element displacement parameter set is generated.

[0051] In a second aspect, a computing device includes:

[0052] One or more processors;

[0053] A storage device for storing one or more programs that, when executed by one or more processors, enable the one or more processors to implement the system.

[0054] Thirdly, a computer-readable storage medium storing a program that, when executed by a processor, implements the system.

[0055] The above-described solution of the present invention has at least the following beneficial effects:

[0056] By employing feature point matching algorithms and spatial transformation parameter optimization techniques, registration errors in sparsely characterized areas such as smooth bare rock and fractured terrain are effectively controlled, resulting in lower errors compared to existing technologies. Through unified spatial benchmark transformation and multi-period data fusion, consistency of monitoring data across time series is ensured, providing a reliable data foundation for displacement analysis.

[0057] A refined land cover classification algorithm was used to accurately separate vegetation cover, artificial structures, and natural surfaces, improving classification accuracy in areas with low to medium vegetation cover and effectively solving the problem of distorted terrain feature extraction in low shrub areas. The optimized terrain feature dataset provides a high-quality data source for building digital terrain surface models, improving model accuracy. Attached Figure Description

[0058] Figure 1 This is a schematic diagram of a slope displacement monitoring data processing system based on UAV lidar provided in an embodiment of the present invention.

[0059] Figure 2 This is a flowchart illustrating the process of extracting and matching stable terrain feature points in the observation dataset, using a feature point matching algorithm to calculate spatial transformation parameters, accurately registering the optimized lidar observation data, completing the conversion and fusion of all data to a unified spatial reference, and generating a slope three-dimensional observation dataset. Detailed Implementation

[0060] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0061] like Figure 1 As shown, an embodiment of the present invention proposes a slope displacement monitoring data processing system based on UAV lidar, comprising:

[0062] The slope displacement monitoring data processing system based on UAV lidar includes:

[0063] The acquisition module is used to acquire the raw lidar observation dataset;

[0064] The processing module is used to preprocess the original lidar observation dataset to obtain an optimized lidar observation dataset.

[0065] The calculation module is used to extract and match stable terrain feature points in the observation dataset, calculate spatial transformation parameters using a feature point matching algorithm, accurately register the optimized lidar observation data, complete the conversion and fusion of all data to a unified spatial reference, and generate a slope three-dimensional observation dataset.

[0066] The generation module is used to establish an directional reference frame based on a slope three-dimensional observation dataset, with two preset stable reference benchmarks, and to generate adaptive correction factors according to the topological relationships and structural characteristics of each region within the frame.

[0067] The classification module optimizes the land cover classification parameters based on an adaptive correction factor to achieve accurate separation of vegetation cover, artificial structures and natural land surface, and generate an optimized terrain feature dataset.

[0068] The building module is used to perform spatial interpolation on the terrain feature dataset to build a digital terrain surface model;

[0069] The analysis module is used to calculate the displacement vector field of the slope surface based on the digital terrain surface model, through spatial similarity analysis and gradient descent algorithm for iterative optimization, identify potential sliding surfaces and abnormal deformation areas, and generate displacement field solution results.

[0070] The output module is used to input the displacement field calculation results into the risk assessment model and perform quantitative evaluation of slope stability through a multi-source data fusion analysis platform.

[0071] In this embodiment of the invention, by employing a feature point matching algorithm and spatial transformation parameter optimization technology, the registration error in sparsely characterized areas such as smooth bare rock and fractured terrain is effectively controlled, which is lower than that of the prior art. By unifying spatial benchmark transformation and fusing data from multiple periods, the consistency of cross-time series monitoring data is ensured, providing a reliable data foundation for displacement analysis.

[0072] By using a fine-grained land cover classification algorithm, the precise separation of vegetation cover, artificial structures and natural land surface is achieved, improving the classification accuracy in areas with low to medium vegetation cover. This effectively solves the problem of distorted extraction of terrain features in low shrub areas. The optimized terrain feature dataset provides a high-quality data source for the construction of digital terrain surface models, improving the accuracy of the models.

[0073] In a preferred embodiment of the present invention, the original lidar observation dataset is preprocessed to obtain an optimized lidar observation dataset, including:

[0074] The original lidar observation dataset is noise filtered, and discrete noise points are removed based on statistical outlier analysis to generate a preliminary clean dataset. The preliminary clean dataset is then subjected to intensity correction, and intensity distortion caused by measurement geometry effects is eliminated through distance and incident angle normalization to generate an intensity-normalized dataset. Specifically, when noise filtering the original lidar observation dataset, for each 3D point in the dataset, all neighboring points with a preset radius of 0.3 meters to 1 meter centered on that 3D point are selected. The Euclidean distance between the 3D point and each neighboring point is calculated. Based on all Euclidean distances, the mean distance and standard deviation are calculated. Neighboring points whose Euclidean distances to each neighboring point are greater than the sum of the mean distance and three times the standard deviation are marked as potential noise points. After removing all potential noise points, the preliminary clean dataset is generated. When performing intensity correction on the initial cleaned dataset, the actual measured distance between the transmitter and each three-dimensional point is obtained through the ranging module built into the lidar sensor. This distance ranges from 50 meters to 200 meters. The incident angle when the laser beam reaches the three-dimensional point is obtained through the angle measurement unit built into the sensor. The incident angle ranges from 15 degrees to 60 degrees.

[0075] When constructing the intensity normalization model, the correlation between the original intensity value and the actual measured distance and incident angle is determined. The original intensity value is used as the dependent variable, and the reciprocal of the square of the actual measured distance and the cosine of the incident angle are used as independent variables to establish a multivariate regression model framework. During the model training phase, a standard reference target with known reflectivity is selected within the slope area. Measurement points are set at 10-meter intervals within a distance range of 50 to 200 meters. At each measurement point, the laser emission angle is adjusted within an incident angle range of 15 to 60 degrees at 5-degree intervals. The original intensity data of the standard target is collected to form a training sample set. Data fitting is performed on the training sample set, and distance correction coefficients and angle correction coefficients are calculated. The distance correction coefficient ranges from [value missing]. The values ​​of the angle correction coefficient range from 0.9 to 1.1, and the coefficients are optimized using the least squares method to keep the deviation between the model output value and the theoretical reflection intensity of the standard target within 5%. During the calculation of the least squares optimization coefficients, an error function is established, and the square of the difference between the model output correction intensity value and the theoretical reflection intensity of the standard target for each training sample is used as a single sample error term. The model output correction intensity value is calculated based on the initial distance correction coefficient and the angle correction coefficient. The single sample error terms of all training samples are summed to obtain the total error function. The expression of the total error function is the sum of the error terms of each sample. The sample error term is equal to the square of the model output correction intensity value and the theoretical reflection intensity of the standard target.

[0076] The partial derivatives of the distance correction coefficient and the angle correction coefficient are calculated for the total error function. Based on the extreme condition that the partial derivatives are zero, a system of two linear equations in two variables is established. The system of equations contains two unknowns: the distance correction coefficient and the angle correction coefficient. The coefficients of the equations are calculated from the actual measured distance, incident angle and original intensity value in the training sample set. The system of two linear equations in two variables is solved by matrix solving method to obtain the initial optimized values ​​of the distance correction coefficient and the angle correction coefficient. Then, the initial optimized value is substituted into the total error function to calculate the current total error. If the total error is greater than the preset error threshold, which is set to 5% of the theoretical reflection intensity of the standard target, the distance correction coefficient and angle correction coefficient are adjusted according to the preset step size (the step size is 0.01) starting from the initial optimized value. The total error is recalculated and the above process of finding partial derivatives and solving the system of equations is repeated until the total error is less than or equal to the error threshold. The distance correction coefficient and angle correction coefficient obtained at this time are the final optimization results. In the model implementation process, for each three-dimensional point in the preliminary cleaned dataset, its actual measured distance and incident angle are extracted and substituted into the intensity normalization model to calculate the corrected intensity value, that is, the corrected intensity value = original intensity value × distance correction coefficient × (square of the actual measured distance reference value / square of the actual measured distance) × angle correction coefficient × cos(incident angle), where the actual measured distance reference value is set to 100 meters. The intensity value adjusted by this model eliminates the distortion caused by the difference in actual measured distance and the change in incident angle, and generates an intensity normalized dataset.

[0077] Coordinate correction is performed on the intensity-normalized dataset. Based on sensor pose parameters and inertial measurement unit (IMU) data, geometric deviations caused by attitude changes of the UAV platform are compensated to generate a geometric correction dataset. Specifically, when performing coordinate correction on the intensity-normalized dataset, the three-dimensional position parameters of the lidar sensor at the time of data acquisition are obtained from the UAV's onboard global navigation satellite system. The position parameters have an accuracy of centimeters. The pitch angle, roll angle, and yaw angle of the sensor at the same time are obtained from the IMU. The attitude parameters have an accuracy of 0.1 degrees. A coordinate transformation matrix is ​​constructed based on the translation transformation of the position parameters and the rotation matrix of the attitude parameters. This matrix is ​​used to perform transformation calculations on the original coordinates of each three-dimensional point in the intensity-normalized dataset to compensate for the geometric deviations caused by attitude changes during UAV flight, thus generating a geometric correction dataset.

[0078] The geometric correction dataset is sampled and optimized, and the point density is adaptively adjusted based on slope topographic features. Redundant data points are removed to generate an optimized lidar observation dataset. Specifically, during the sampling and optimization of the geometric correction dataset, the 3D point cloud data in the geometric correction dataset is partitioned into 10m × 10m grids. The slope value of the 3D points in each partition is calculated based on the elevation difference and horizontal distance between adjacent points. The curvature value of each partition is also calculated. The target point density for different partitions is set according to the slope and curvature values, where the density is set at a ratio of slope greater than 30 degrees or curvature greater than 0.05. The target point density is set to 20 to 30 points / m² for zones with a slope of 15 to 30 degrees and a curvature of 0.02 to 0.05 meters, and 10 to 20 points / m² for zones with a slope of less than 15 degrees and a curvature of less than 0.02 meters. The target point density is set to 5 to 10 points / m² for zones with a slope of less than 15 degrees and a curvature of less than 0.02 meters. The 3D points in each zone are uniformly sampled and filtered according to the target point density. Key feature points are retained and redundant data points that exceed the target density are removed to generate an optimized lidar observation dataset.

[0079] In this embodiment of the invention, discrete noise is identified through statistical characteristics, with a value range adaptable to different slope point cloud densities, effectively preserving terrain feature points while reducing noise interference. This step corrects intensity distortion through physical characteristic modeling, with a value range covering the conventional monitoring altitude and angle of UAVs, which can improve the consistency of ground feature intensity characteristics in different areas.

[0080] like Figure 2 As shown, in another preferred embodiment of the present invention, stable terrain feature points in the observation dataset are extracted and matched, spatial transformation parameters are calculated using a feature point matching algorithm, and the optimized lidar observation data is precisely registered to complete the conversion and fusion of all data to a unified spatial reference, generating a slope three-dimensional observation dataset, including:

[0081] Stable terrain feature points are extracted from the observation data. Based on the curvature and normal vector change rate indices, surface protrusions, rock strata boundaries, and abrupt topographic changes are selected as high-reliability feature points to generate an initial feature point set. Specifically, when extracting stable terrain feature points from the observation data, for each 3D point in the optimized lidar observation dataset, a neighborhood point set with a radius of 0.5 to 2 meters centered on that point is selected. A local quadratic surface is fitted using the neighborhood point set, and the curvature value of that point is calculated. The curvature value is calculated using the flattening of the principal curvatures of the surface. The mean value is calculated, ranging from 0 to 0.5 per meter. The angle between the normal vector of the point and the normal vector of all points in the neighborhood is calculated. The rate of change of the normal vector is obtained by the ratio of the angle difference to the distance between the points. The rate of change of the normal vector ranges from 0 to 30 degrees / meter. The curvature threshold is set to 0.1 per meter and the rate of change of the normal vector threshold is set to 15 degrees / meter. Points with curvature values ​​greater than or equal to the curvature threshold or the rate of change of the normal vector greater than or equal to the rate of change of the normal vector threshold are selected and marked as surface protrusions, rock layer boundaries, or abrupt topographic changes, thus generating an initial feature point set.

[0082] When using the least squares method to fit the surface to the sample points, a quadratic surface model to be fitted is determined. This model includes quadratic terms for the x-coordinate, quadratic terms for the y-coordinate, the intersection of x and y, linear terms for the x-coordinate, linear terms for the y-coordinate, and a constant term. This model describes the terrain surface morphology of the local area where the sample points are located. When selecting sample points, a neighborhood with a radius of 0.5 to 2 meters is defined, centered on the target 3D point. All 3D points within this neighborhood are extracted as fitting sample points. The number of sample points is controlled between 50 and 200 to ensure sufficient coverage of terrain feature information. When constructing the error function, the squared difference between the actual z-coordinate value of each sample point and the predicted z-coordinate value calculated through the surface model is used as a single sample error term. The single sample error terms of all sample points are summed to obtain the total error function, which reflects the degree of deviation between the model's predicted value and the actual terrain.

[0083] Calculate the partial derivatives of each coefficient in the model for the total error function. Establish a system of linear equations based on the extreme condition that the partial derivatives are zero. The unknowns in the system of equations are the coefficients of the surface model. The coefficient matrix of the equations consists of the x-coordinates, y-coordinates, and constant terms of the sample points. The constant term vector consists of the sum of the products of the actual z-coordinates and the combined coordinate values ​​of the sample points. Solve the system of linear equations by matrix inversion to obtain the initial values ​​of the coefficients of the surface model. Substitute the initial coefficients into the surface model to calculate the predicted z-coordinates of each sample point. Calculate the residuals between the actual z-coordinates and the predicted z-coordinates. Calculate the standard deviation based on the residuals. If the standard deviation is greater than a preset threshold (set to 0.05 m), remove outlier sample points whose absolute residual value is greater than twice the standard deviation. Reselect the remaining sample points and repeat the above fitting process until the residual standard deviation is less than or equal to the threshold or the number of iterations reaches a preset upper limit (set to 5 iterations). The coefficients obtained at this point are the final coefficients of the final surface model.

[0084] The total error function is constructed using the least squares method. The total error function E is the sum of squared errors between the actual and predicted z-coordinates of all sample points, and its calculation formula is as follows: Where n is the total number of sample points, Let z be the actual z-coordinate value of the i-th sample point. Let be the predicted coordinates of the i-th sample point; Substituting into the above equation, the total error function can be expanded to E= ,in, , Let i be the planar coordinates of the i-th sample point. , , , , , Let be the coefficients to be determined for the quadratic polynomial model; after expansion, it can be further expressed as E= + ;in, It is the linear contribution term of the x-coordinate of the i-th sample point; It is the linear contribution term of the ordinate of the i-th sample point; It is the nonlinear contribution term of the x-coordinate of the i-th sample point; It is the coupling contribution term of the horizontal and vertical coordinates of the i-th sample point; It is the nonlinear contribution term of the ordinate of the i-th sample point.

[0085] , These represent the x-coordinate and y-coordinate (horizontal coordinates) of the i-th sample point in three-dimensional space, respectively, which are used to describe the planar position of the sample point. This represents the actual z-coordinate (elevation value) of the i-th sample point in three-dimensional space, i.e., the true terrain height of the sample point; This represents the predicted z-coordinate (predicted elevation) of the i-th sample point calculated using the quadratic surface model, used to compare with the actual z-coordinate. The comparison was used to evaluate the model's fit. The constant term coefficients in the quadratic surface model represent one of the fundamental parameters for model fitting. In the quadratic surface model, the first-order term of the x-coordinate is represented as ( The coefficient of ) is used to characterize the linear variation trend of terrain in the x-direction; This represents the first-order term of the y-coordinate in the quadratic surface model. The coefficient of ) is used to characterize the linear variation trend of terrain in the y direction; This represents the quadratic term of the x-coordinate in the quadratic surface model. The coefficient of ) is used to characterize the nonlinear curvature of the terrain in the x-direction; This represents the x-y intersection term in the quadratic surface model. • The coefficients of ) are used to characterize the nonlinear features of terrain under the interaction of the x and y directions; This represents the quadratic term of the y-coordinate in the quadratic surface model. The coefficient of ) is used to characterize the nonlinear curvature of the terrain in the y-direction.

[0086] To address the need for fitting local topographic data from lidar observations in slope monitoring, a quadratic surface model is constructed as a mathematical carrier to describe local topographic morphology. This model must fully reflect the nonlinear relationship between the x, y, and z coordinates to adapt to both rock and soil areas of the slope. The model structure is set as a complete quadratic surface containing six terms, specifically covering the quadratic terms of the x and y coordinates, the x-y intersection term, the x-y linear term, the y-y linear term, and a constant term. Each term corresponds to an independent coefficient. The x² term coefficient is used to characterize the degree of topographic curvature along the x direction, the y² term coefficient is used to characterize the degree of topographic curvature along the y direction, the xy intersection term coefficient is used to characterize the degree of topographic distortion in the x and y planes, the x-y linear term coefficient and the y-y linear term coefficient represent the topographic tilt trend along the x and y directions, respectively, and the constant term represents the elevation datum of the model at the origin.

[0087] The selection of sample points needs to be determined based on the complexity of the slope terrain. A circular neighborhood with a radius ranging from 0.5 meters to 2 meters is defined, centered on the target 3D point. In rocky areas, where the terrain is relatively regular, the neighborhood radius is 0.5 meters to 1.2 meters. In soily areas, where the loose soil contains boulders and the terrain is more complex, the neighborhood radius is 1.2 meters to 2 meters. This ensures that the neighborhood contains a sufficient number of sample points to reflect the local terrain features. The number of sample points is controlled between 50 and 200. For rocky areas, the number of sample points per neighborhood is 50 to 120, and for soily areas, it is 120 to 200. If the number of sample points in a neighborhood is less than 50, the neighborhood radius is expanded until the required number of sample points is reached. If the number of sample points exceeds 200, a uniform sampling method is used to remove some redundant points (retaining one point every two points) to avoid excessive data increasing the computational load while ensuring the representativeness of the sample points.

[0088] The selected sample points are arranged in coordinate format as (x1, y1, z1), (x2, y2, z2)...(x... n y n , z n The structured data consists of x and y coordinates, which are two-dimensional position coordinates in the horizontal plane, and z coordinate, which is the actual elevation value measured by lidar. The coordinates of all sample points are converted to the geodetic coordinate system uniformly used for slope monitoring to ensure the consistency of coordinate reference and lay the foundation for subsequent error calculation and coefficient solution.

[0089] Feature descriptors are constructed from the initial feature point set, generating feature description vectors based on local surface geometric attributes to obtain a standardized feature description set. Similarity matching is then performed on the standardized feature description set, and mismatched point pairs are eliminated through a bidirectional consistency check to generate a reliable feature point correspondence set. Specifically, when constructing feature descriptors from the initial feature point set, for each feature point, a local neighborhood with a radius of 1 to 3 meters is selected centered on that point. The average curvature, normal vector direction, elevation standard deviation, and intensity mean of all points within the neighborhood are calculated, and the geometric attributes are arranged in a preset order to form a basic feature vector. The basic feature vector is then standardized, converting each attribute value to the range of 0 to 1. Linear scaling is used to eliminate dimensional differences between different attributes, forming a standardized feature description vector of length 128 dimensions. The description vectors of all feature points constitute the standardized feature description set.

[0090] Based on the reliable feature point correspondence set, the final spatial transformation parameters are calculated, and a coordinate transformation parameter set is generated. The coordinate transformation parameter set is then used to perform a rigid transformation on all lidar observation data, uniformly converting them to the reference coordinate system, and generating a dataset with a unified spatial reference. Specifically, when calculating the final spatial transformation parameters based on the reliable feature point correspondence set, point pairs in the reliable feature point correspondence set are used as input to construct a spatial transformation model composed of a rotation matrix and a translation vector. The initial transformation parameters are calculated using the three-point method. An iterative optimization method is used to minimize the distance error of corresponding points after transformation. In each iteration, the transformed Euclidean distance of all corresponding points is calculated. Point pairs whose distance is greater than the sum of the average distance and twice the standard deviation are marked as outliers and temporarily removed. The transformation parameters are recalculated based on the remaining inliers. The iterative process is repeated until the change in transformation parameters between two consecutive iterations is less than 0.01 radians and 0.05 meters. The resulting rotation matrix and translation vector form the coordinate transformation parameter set. When performing rigid transformation on all periods of lidar observation data using the coordinate transformation parameter set, the three-dimensional point coordinates in each period's data are substituted into the rotation matrix and translation vector in the coordinate transformation parameter set. The angular deviation between different periods of data is eliminated through rotation operations, and the spatial position benchmark of the data is unified through translation operations. The benchmark coordinate system is selected from the coordinate system of the first period of observation data. The coordinates of all data points after transformation are converted to this benchmark coordinate system, generating a dataset with a unified spatial benchmark.

[0091] The construction, training, and implementation process of the spatial transformation model: When constructing the spatial transformation model, a three-dimensional rigid transformation model is used as the basic framework. This model contains two core parameters: rotation and translation. The rotation parameter describes the spatial rotation relationship through three Euler angles (pitch angle, roll angle, and yaw angle), while the translation parameter describes the spatial position adjustment relationship through three-dimensional coordinate offsets (Δx, Δy, Δz). The mathematical expression of the model is that the three-dimensional coordinates in the target coordinate system are equal to the original coordinates after transformation by the rotation matrix and the superposition of the translation vector, thereby achieving spatial position alignment of LiDAR data from different periods.

[0092] During the model training phase, a reliable set of feature point correspondences is used as the input sample. 30 to 100 pairs of feature points are randomly selected from the correspondence set as training sample pairs. The sample pairs should be evenly distributed in the rocky area, soil area and boundary area of ​​the slope to avoid concentration in a single terrain unit. When establishing the error function, the square of the Euclidean distance between the predicted coordinates of each pair of feature points calculated by the transformation model and the actual target coordinates is used as a single sample error term. The error terms of all sample pairs are summed to obtain the total error function, which reflects the overall fitting effect of the transformation model. The partial derivatives of the rotation and translation parameters of the total error function are calculated separately. A nonlinear equation system is established based on the extremum condition that the partial derivatives are zero. The equation system is solved using the Gauss-Newton iteration method: the initial value of the rotation angle is set to 0 degrees, and the initial value of the translation vector is set to 0. The parameter corrections are calculated by substituting these values ​​into the equation system. After updating the rotation matrix and translation vector, the total error is recalculated. This iterative process is repeated until the difference in total error between two iterations is less than 0.01 meters or the number of iterations reaches 20, yielding the preliminary optimized transformation parameters. Outlier cleanup is performed on the preliminary optimized parameters. The Euclidean distance of the residuals after transformation is calculated for all feature point pairs. The mean and standard deviation of the residuals are statistically analyzed. Feature point pairs with residuals greater than the mean plus twice the standard deviation are marked as outlier pairs. The remaining feature point pairs are then discarded, and the above training process is repeated to generate the final optimized spatial transformation parameter set. The rotation angle accuracy is controlled within 0.01 degrees, and the translation vector accuracy is controlled within 0.05 meters. During model implementation, the optimized rotation matrix and translation vector are applied to all periods of LiDAR observation data. The original coordinates of each 3D point are rotated and translated to the reference coordinate system. After the transformation, 50 to 100 feature points from different areas of the slope are randomly selected to verify the registration accuracy. The average distance error of the registered corresponding points is calculated to ensure that the error is less than 0.1 meters. If the error exceeds the error threshold, training sample pairs are reselected for model training and optimization until the accuracy requirements are met.

[0093] The process involves fusing datasets with a unified spatial reference, eliminating redundant points in overlapping areas, and generating a complete three-dimensional slope observation dataset. Specifically, this includes: dividing the dataset into 5m x 5m grids, counting the number of 3D points within each grid, and when the number of points in a grid exceeds a preset density threshold (set to 50 points per grid), a uniform sampling method is used to retain the point cloud, i.e., one point is retained every 0.5 meters according to the spatial distribution of points within the grid, while redundant points exceeding the density threshold are removed; when the number of points in a grid is below the density threshold, all points are retained to ensure terrain integrity. The fused data covers the entire slope area without overlap or redundancy, generating a complete three-dimensional slope observation dataset.

[0094] In this embodiment of the invention, a standardized feature description set is constructed by using local surface geometric properties, and mismatched point pairs are eliminated by combining a two-way consistency check, which can reduce matching errors caused by similar features in complex terrain. Stable terrain feature points are screened from the observation data based on curvature and normal vector change rate indices, which can accurately identify geometrically distinctive feature points such as surface protrusions, rock layer boundaries, and abrupt terrain changes, while eliminating unstable points susceptible to environmental interference.

[0095] In a preferred embodiment of the present invention, based on a slope three-dimensional observation dataset, two stable reference benchmarks are preset, an oriented reference frame is established, and an adaptive correction factor is generated according to the topological relationships and structural characteristics of each region within the frame, including:

[0096] From the slope 3D observation dataset, a first spatial reference point located in the stable bedrock region at the top of the slope and a second spatial reference point located in the stable bedrock region at the toe of the slope are selected to generate an initial reference axis, specifically including:

[0097] Using the initial reference axis as the reference direction, a regional coordinate system is constructed with the first spatial reference point as the origin and the direction of the reference line as the main axis, generating an directional reference frame. Specifically, this includes: conducting a comprehensive preliminary geological survey of the slope of a specific section of the highway in the mountainous area, clarifying the bedrock distribution in the top and toe areas of the slope, determining the area at the top of the slope without the risk of bedding-parallel sliding and with bedrock integrity meeting the requirements as the candidate area for the first spatial reference point, and determining the area at the toe of the slope without the accumulation of loose silty clay and with no risk of bedrock weathering and erosion as the candidate area for the second spatial reference point.

[0098] A three-dimensional observation dataset for the slope was acquired using a combination of UAV oblique photogrammetry and terrestrial 3D laser scanning. UAV oblique photogrammetry was used to obtain overall 3D topographic data of the slope, while terrestrial 3D laser scanning was used to obtain high-precision point cloud data for candidate areas at the top and toe of the slope. The point cloud data of the candidate areas was then denoised to remove invalid point cloud data caused by tree obstruction and dust interference, retaining only valid point cloud data reflecting the actual morphology of the bedrock. From the processed valid point cloud data of the candidate area at the top of the slope, three non-collinear points were selected. The coordinates of these three points were measured in the field using a total station, and the average coordinate of these three points was calculated as the 3D coordinates of the first spatial reference point. This ensures that the first spatial reference point is located in the stable bedrock area at the top of the slope and can represent the spatial location benchmark of that area. Simultaneously, three non-collinear points are selected from the effective point cloud data of the candidate area at the slope toe after processing. The coordinates of these three points are measured in the field using a total station, and the average coordinate of these three points is calculated as the three-dimensional coordinate of the second spatial reference point. This ensures that the second spatial reference point is located in the stable bedrock area at the slope toe and can represent the spatial location benchmark of the area. The three-dimensional coordinates of the first spatial reference point and the second spatial reference point are connected by a line using three-dimensional coordinate calculation software to form a straight line connecting the two reference points. This straight line is the initial reference axis.

[0099] Based on the directional reference frame, and according to the geological structure characteristics and topographic change trends of the slope, the monitoring area is divided into several independent analysis units, generating unit division results. Specifically, the first spatial reference point is set as the origin of the regional coordinate system, and the direction of the initial reference axis from the first spatial reference point to the second spatial reference point is set as the principal axis direction of the regional coordinate system. This principal axis direction is consistent with the length direction of the slope to cover the entire slope monitoring area. The magnetic north direction of the area where the slope is located is measured by a geological compass, and the horizontal coordinate axis direction perpendicular to the principal axis direction is determined by combining the topographic slope of the slope. This horizontal coordinate axis direction is parallel to the transverse section direction of the slope and is used to reflect the spatial position of the slope width direction.

[0100] The vertical axis direction of the regional coordinate system is determined, and this vertical axis direction is perpendicular to the local horizontal plane. The reference height of the local horizontal plane is obtained by measuring with a level instrument to ensure that the vertical axis direction accurately reflects the height direction of the slope. The three coordinate axes of the regional coordinate system are calibrated, and the correspondence between the unit length of each coordinate axis and the actual physical length is clarified. The unit length is set to a length unit that can accurately reflect the small deformations of the slope to meet the monitoring accuracy requirements. Then, all three-dimensional data in the slope stereoscopic observation dataset are transformed to the regional coordinate system. The coordinate values ​​of the original observation data are converted to coordinate values ​​in the regional coordinate system through a coordinate transformation algorithm to ensure that all monitoring data are processed in the same coordinate system. The origin position, coordinate axis direction, unit length, and data transformation rules of the regional coordinate system are organized to form an orientation reference frame. This orientation reference frame can accurately correspond to the actual spatial location and topographic features of the slope, providing a unified spatial reference for subsequent analysis unit division.

[0101] Geological structural feature data of the slope was extracted within the directional reference frame. This data included the distribution range of the upper moderately weathered sandstone interbedded with mudstone strata, the angle between the strata dip and the slope aspect, the distribution range of the lower Quaternary residual colluvial silty clay, the degree of soil looseness, and the distribution of boulders. At the same time, topographic change trend data of the slope was extracted, including the slope change of different areas of the slope, the location distribution of the slope shoulder, slope waist, and slope toe, and the topographic undulation characteristics. Secondly, based on the extracted geological structural feature data, the demarcation boundary was determined. The boundary between the upper moderately weathered sandstone interbedded with mudstone strata and the lower silty clay strata was set as the longitudinal demarcation boundary, and the boundary between the area with potential for bedding-parallel sliding and the area without potential for bedding-parallel sliding was set as the lateral demarcation boundary.

[0102] The division intervals were determined by combining topographic change trend data. In silty clay areas with steep slopes and prone to landslides, and in sandstone-mudstone areas with prominent bedding-parallel sliding risks, smaller division intervals were used to increase the density of unit divisions. In areas with gentle slopes and stable geological conditions, larger division intervals were used to reduce redundant units. According to the determined division boundaries and intervals, the entire slope monitoring area was divided into multiple independent analysis units within the directional reference frame. The spatial range of each independent analysis unit could completely encompass a specific geological structure feature and topographic change feature, ensuring the consistency of geological and topographic conditions within each unit. Each independent analysis unit was numbered, and the coordinates of the four vertices of each unit within the directional reference frame, as well as the corresponding geological structure type and topographic feature description, were recorded. The information was then compiled to form the unit division results.

[0103] The topographic structure attributes of each analysis unit are extracted to generate a set of unit topographic feature descriptions. Based on the relative spatial position relationship between each analysis unit and two spatial reference points, the azimuth weight and distance weight of each unit relative to the reference line are calculated to generate a set of unit weight distributions. Specifically, for each independent analysis unit in the unit division results, the topographic structure attributes of the unit are extracted using three-dimensional topographic analysis software. The extracted topographic structure attributes include the unit's average slope, maximum slope, minimum slope, average aspect, range of aspect variation, surface relief, and the corresponding stratum type. If the unit belongs to the lower silty clay stratum, the density grade and proportion of boulders in the silty clay within the unit also need to be extracted. These attribute data are obtained by statistical analysis of the point cloud data of the corresponding units in the slope three-dimensional observation dataset. Then, the topographic structure attributes of each independent analysis unit are organized according to a preset format. Each unit corresponds to one attribute record, which includes the unit number, average slope, maximum slope, minimum slope, average aspect, range of aspect variation, surface relief, stratum type, density level, and proportion of block content. The attribute records of all units are summarized to form a unit topographic feature description set.

[0104] To calculate the azimuth weight of each independent analysis unit relative to two spatial reference points, first, the coordinates of the geometric center point of each unit are determined within the orientation reference frame. The angle between this geometric center point and the initial reference axis is calculated, which is the azimuth angle of the unit. An azimuth weight mapping rule is set according to the slope's risk characteristics. Units located within the azimuth angle range of the bedding slip hazard section and the azimuth angle range of the silty clay prone-to-collapse section are assigned higher azimuth weight values, while units located within the azimuth angle range of the geologically stable area are assigned lower azimuth weight values. The azimuth weight of each unit is calculated according to this mapping rule. The distance weights of each independent analysis unit relative to the two spatial reference points are calculated separately. The straight-line distance from the geometric center point of each unit to the first spatial reference point and the straight-line distance to the second spatial reference point are calculated separately. The distance weight calculation rules are set according to the correlation between distance and reference point stability. Units closer to the two spatial reference points are assigned higher distance weight values ​​because they are more affected by the stability of the reference points, while units farther away from the two spatial reference points are assigned lower distance weight values. The distance weight of each unit is calculated separately according to the calculation rules. The unit number, orientation weight value and distance weight value of each independent analysis unit are organized to form a unit weight distribution set.

[0105] The method integrates the unit topographic feature description set and the unit weight distribution set to calculate the comprehensive topographic stability index of each analysis unit and generate a unit stability evaluation set. Specifically, this includes: quantifying the various topographic structural attributes in the unit topographic feature description set; setting quantitative scoring standards based on the influence of each attribute on slope stability, where a larger average slope, a smaller average slope aspect and rock stratum dip angle, a larger surface undulation, a lower silty clay density grade, a higher proportion of boulder, and a higher score for moderately weathered sandstone interbedded with mudstone strata than silty clay strata; scoring each topographic structural attribute of each independent analysis unit according to this quantitative scoring standard to obtain the attribute score value for each unit; extracting the azimuth weight value and distance weight value of each independent analysis unit from the unit weight distribution set; and merging the azimuth weight value and distance weight value according to a preset ratio to obtain the comprehensive weight value of each unit. Considering that the risk of slope bedding slip and collapse is mainly related to azimuth, the azimuth weight value has a higher proportion in the comprehensive weight value calculation than the distance weight value.

[0106] Each attribute score of each independent analysis unit is multiplied by the overall weight value of that unit to obtain the weighted score for each attribute. The weighted scores of all attributes for each unit are then summed to obtain the overall terrain stability index for that unit. The larger the value of the overall index, the better the terrain stability of the unit, and the smaller the value, the worse the terrain stability of the unit. The unit number and the corresponding overall terrain stability index of each independent analysis unit are arranged in order of unit number to form a unit stability evaluation set.

[0107] Based on the unit stability evaluation set, stability indices are converted into adjustment coefficients for land feature classification parameters, generating adaptive correction factors. Specifically, this includes: clarifying the types of land feature classification parameters involved in slope monitoring, including slope surface displacement monitoring parameters, deep slope displacement monitoring parameters, soil pore water pressure monitoring parameters, and surface crack development degree monitoring parameters. These land feature classification parameters are key monitoring parameters reflecting slope stability, and their monitoring accuracy directly affects the slope instability risk assessment results; and setting adjustment coefficient values ​​based on the comprehensive topographic stability index values ​​of each independent analysis unit in the unit stability evaluation set. For units with a stable stability level, since their geological and topographic stability monitoring parameters are less affected by environmental interference, their corresponding land feature classification parameter adjustment coefficients are set to smaller values ​​to reduce unnecessary correction operations.

[0108] For units with a relatively stable stability level, due to the slight risk of stability fluctuations, the adjustment coefficients for their corresponding land cover classification parameters are set to medium values ​​to achieve moderate correction. For units with an unstable stability level, due to their susceptibility to bedding slips or landslides and the susceptibility of monitoring parameters to complex environmental influences, the adjustment coefficients for their corresponding land cover classification parameters are set to larger values ​​to ensure the accuracy of the monitoring parameters. For each independent analysis unit, based on the stability level corresponding to its comprehensive topographic stability index and the preset adjustment coefficient value rules, the adjustment coefficients for each type of land cover classification parameter are calculated to ensure that each land cover classification parameter has an adjustment coefficient that matches the unit's stability. The unit number, land cover classification parameter names, and corresponding adjustment coefficients for each independent analysis unit are compiled and summarized to form an adaptive correction factor.

[0109] In this embodiment of the invention, by using a dual-benchmark orientation and zonal weighting mechanism, the unit size and weight range adapt to the differences in composite slope terrain, which can improve the correction accuracy in complex terrain areas. By adjusting the parameters related to stability, it provides an adaptive parameter optimization basis for subsequent land feature classification processing, reducing classification bias caused by terrain variations.

[0110] In a preferred embodiment of the present invention, the land cover classification processing parameters are optimized according to an adaptive correction factor to achieve accurate separation of vegetation cover, artificial structures, and natural land surfaces, generating an optimized terrain feature dataset, including:

[0111] Based on the adaptive correction factor, the parameter adjustment coefficients of each analysis unit are analyzed to generate a set of unit parameter adjustments. Based on this set, the classification threshold parameters in the land cover classification algorithm are adaptively adjusted according to the terrain stability characteristics of different analysis units to generate an optimized classification parameter set. Specifically, when optimizing the land cover classification processing parameters based on the adaptive correction factor, the parameter adjustment coefficients of each analysis unit need to be matched one by one according to the analysis unit. Each analysis unit corresponds to a unique parameter adjustment coefficient, with a value range of 0.6 to 1.4. Units with high terrain stability have coefficients greater than 1.0, while units with low terrain stability have coefficients greater than 1.0. If the element correspondence coefficient is less than 1.0, a unit parameter adjustment set containing all unit adjustment coefficients is generated. When adjusting the classification threshold parameters based on the unit parameter adjustment set, the initial reflectance intensity threshold range is set to 20 to 60 for vegetation cover classification, 60 to 100 for artificial structure classification, and 10 to 30 for natural surface classification. For high-stability units, the threshold interval between vegetation and natural surface is increased by 10% to 20% through parameter adjustment coefficients, and for low-stability units, the threshold interval between artificial structure and natural surface is reduced by 5% to 15%, generating an optimized classification parameter set for each analysis unit.

[0112] Preliminary land cover classification is performed on the slope 3D observation dataset based on the optimized classification parameter set to generate initial classification results. The initial classification results are then processed to generate optimized classification results. Specifically, during the preliminary land cover classification based on the optimized classification parameter set, the reflection intensity value and local slope value of each 3D point in the slope 3D observation dataset are extracted point by point. Points with reflection intensity values ​​falling within the vegetation threshold range and slope less than 25 degrees are marked as vegetation-covered points; points with reflection intensity values ​​falling within the artificial structure threshold range and slope less than 15 degrees are marked as artificial structure points; the remaining points are marked as natural surface points, generating initial classification results. During the processing of the initial classification results, the number of 3D points in each classification region is counted. Isolated point groups with fewer than 50 points are identified as classification errors and reclassified to adjacent dominant classification regions. Buffer analysis is performed on the edge points of vegetation-covered areas and artificial structure areas, with a buffer distance set to 0.5 meters. Blurred points within the buffer zone are reclassified according to the majority principle, generating optimized classification results.

[0113] Based on the optimized classification results, spatial distribution information of vegetation-covered areas and artificial structure areas is extracted to generate a ground feature mask dataset. The slope stereoscopic observation dataset is then processed using this dataset to remove data points corresponding to vegetation and artificial structures, generating a three-dimensional point set of exposed ground. Specifically, when generating the ground feature mask dataset based on the optimized classification results, the three-dimensional points of vegetation-covered areas and artificial structure areas are spatially rasterized. The raster size is set to 0.5m × 0.5m, and the rasters containing vegetation or artificial structure points are marked as mask areas, generating a ground feature mask dataset composed of raster coordinates and mask labels. When processing the slope stereoscopic observation dataset using the ground feature mask dataset, the mask labels of the rasters containing the three-dimensional points are checked point by point. All three-dimensional points located within the mask areas are removed, and the three-dimensional points not covered by the mask are retained, generating a three-dimensional point set of exposed ground.

[0114] The three-dimensional point set of exposed land surface is processed to generate an optimized terrain feature dataset. Specifically, when processing the three-dimensional point set of exposed land surface, a neighborhood analysis with a radius of 0.3 meters is used to remove discrete noise points. The criteria for noise points are isolated points with fewer than 5 points in their neighborhood. Then, the point density is calculated by dividing the area into 10-meter × 10-meter grids. For areas with a density lower than 3 points / square meter, interpolation is performed to add points. For areas with a density higher than 20 points / square meter, uniform downsampling is performed. Finally, an optimized terrain feature dataset with a point density controlled between 5 and 15 points / square meter is generated.

[0115] In this embodiment of the invention, the accurate separation of features under complex terrain is achieved by adaptive parameter adjustment. The range of values ​​is adapted to the distribution characteristics of slope vegetation and artificial structures, thereby improving the accuracy of feature classification and effectively eliminating the interference of vegetation and artificial structures on terrain deformation analysis.

[0116] In a preferred embodiment of the present invention, a digital terrain surface model is constructed by spatial interpolation of the terrain feature dataset, including:

[0117] Topographic feature analysis is performed on the optimized topographic feature dataset. Topographic feature points are identified based on local curvature changes and elevation variation coefficients. Ridge lines, valley lines, and abrupt slope changes are selected as key topographic points to generate a candidate topographic feature point set. Density optimization sampling is performed on the candidate topographic feature point set. The sampling density is adaptively adjusted based on the topographic complexity, with sparse sampling in flat areas and dense sampling in complex areas to generate a structured topographic point set. Specifically, when performing topographic feature analysis on the optimized topographic feature dataset, for each 3D point in the dataset, a neighborhood region with a radius of 0.5 meters to 2 meters centered on that point is selected, and 30 to 100 3D points within the neighborhood are extracted as analysis samples. The local curvature value of the point is calculated by surface fitting, with the curvature value ranged from 0 to 0.5 per meter. The ratio of the elevation standard deviation of all points in the neighborhood to the average elevation is calculated as the elevation variation coefficient, with the variation coefficient ranging from 0 to 0.5. Points with a local curvature value greater than 0.1 per meter and an elevation variation coefficient greater than 0.2 are marked as potential terrain feature points. Ridge lines, valley lines, and slope abrupt change points are identified by combining the slope terrain orientation. In the ridge line area, the elevation of more than 3 consecutive potential feature points must be higher than that of the points on both sides. In the valley line area, the elevation of more than 3 consecutive potential feature points must be lower than that of the points on both sides. Slope abrupt change points must have a slope difference of more than 15 degrees between adjacent points. After screening, a candidate terrain feature point set is generated.

[0118] When performing density optimization sampling on the candidate terrain feature point set, the slope area is divided into several sub-regions using a 10m × 10m grid. The average slope and curvature values ​​of each sub-region are calculated. When the average slope of a sub-region is less than 15 degrees and the average curvature is less than 0.05 per meter, it is determined to be a flat terrain area, and the sampling density is set to 5 points / square meter to 10 points / square meter. When the average slope is between 15 degrees and 30 degrees and the average curvature is between 0.05 per meter and 0.1 per meter, it is determined to be a moderately complex terrain area, and the sampling density is set to 10 points / square meter to 20 points / square meter. When the average slope is greater than 30 degrees and the average curvature is greater than 0.1 per meter, it is determined to be a complex terrain area, and the sampling density is set to 20 points / square meter to 30 points / square meter. A uniform sampling method is used to select feature points in each sub-region to ensure that the sampling points uniformly cover the area and that key terrain nodes are not missed, thus generating a structured terrain point set.

[0119] A spatial interpolation function is constructed based on a structured terrain point set to obtain the final interpolation parameters and generate an optimized interpolation function. The optimized interpolation function is then used to calculate the numerical values ​​of regularly distributed elevation points, generating a preliminary terrain surface model. Specifically, this includes: when constructing the spatial interpolation function based on the structured terrain point set, an inverse distance weighted interpolation model is selected as the basic framework, with the spatial location and elevation values ​​of the sampling points as input parameters; during the construction process, the interpolation search radius is set to 5 to 10 meters to ensure that each interpolation calculation point is surrounded by 8 to 20 structured terrain points; and different weights are selected when determining the weight coefficients. Trial calculations were performed using values ​​(ranging from 1 to 5). For each trial calculation, 20% of the structured terrain points were selected as verification points, and the average absolute error between the predicted and actual elevations of the verification points was calculated. When the weighting coefficient was between 2 and 3, the average absolute error was controlled within 0.1 meters. At this point, the weighting coefficient was determined as the final parameter, and an optimized interpolation function was generated. In implementation, the slope area was divided into regular grids with side lengths of 1 to 2 meters. The optimized interpolation function was applied to each grid node, and the node elevation value was obtained by weighting the elevations of the structured terrain points within the search radius, thus generating a preliminary terrain surface model.

[0120] When constructing the optimized interpolation function, an inverse distance weighted interpolation model is used as the basic framework. This model predicts elevation by calculating the distance weights between the interpolation point and surrounding known terrain points; the closer the known points are, the greater their influence on the interpolation results. The slope monitoring area is divided into regular calculation units using a 5m × 5m grid, with the center point of each unit serving as the interpolation point, forming a regularly distributed elevation calculation grid. During the model training phase, three-dimensional points from a structured terrain point set are used as known sample points. A dynamic search radius is set for each interpolation point, adaptively adjusted according to the terrain complexity of the area. In flat terrain, the search radius is set to 15 to 20 meters, while in complex terrain, it is set to 5 to 10 meters, ensuring that 20 to 50 known sample points can be searched around each interpolation point.

[0121] The searched sample points are sorted according to their Euclidean distance from the point to be interpolated. A weight value is calculated for each sample point; the closer the point, the larger the weight value. The weight value decreases non-linearly with increasing distance. The weight decay coefficient is set to 1.5 to 3.0 based on the terrain roughness, with a larger decay coefficient selected for areas with high roughness. When establishing the error evaluation mechanism, 10% of the points are randomly selected from the structured terrain point set as validation samples. These samples do not participate in the interpolation calculation. The elevation values ​​of these validation samples are predicted using an interpolation function. The absolute error between the predicted and actual elevations is calculated, and the average and maximum errors of all validation samples are statistically analyzed. If the average error is greater than 0.1 meters or the maximum error is greater than 0.1 meters, the error is considered a valid error. If the error exceeds 0.3 meters, the search radius and weight decay coefficient are adjusted and recalculated. The search radius is increased or decreased in 2-meter increments, and the decay coefficient is adjusted in 0.5-meter increments. This training process is repeated until the average error is less than or equal to 0.1 meters and the maximum error is less than or equal to 0.3 meters. The resulting search radius range and weight decay coefficient are the final interpolation parameters. During model implementation, the final interpolation parameters are applied to all interpolation points. For the center point of each computational unit, surrounding sample points are selected according to the optimized search radius. The weight values ​​of each sample point are calculated based on the weight decay coefficient. The elevation value of the interpolation point is obtained by weighted averaging, generating a preliminary terrain surface model. During implementation, a local terrain continuity check is performed on each interpolation result. If the elevation difference between adjacent interpolation points exceeds the reasonable difference corresponding to the average slope of the area, the weight parameters are locally fine-tuned to ensure a smooth transition of the terrain surface.

[0122] The initial terrain surface model is smoothed to eliminate local interpolation anomalies, generating an optimized terrain surface model. The optimized model then undergoes integrity verification, filling in data gaps based on terrain continuity principles to generate the final digital terrain surface model. Specifically, during the smoothing of the initial model, a moving window of 3×3 to 5×5 is used to calculate the moving average of the model's elevation data. The elevation value of the center pixel within the window is replaced with the average elevation of all pixels within that window. Pixels with an elevation value greater than 0.2 meters from the average are marked as outliers and replaced with the average of the non-outliers within the window, eliminating abnormal bumps or depressions caused by local interpolation fluctuations. During the integrity verification of the optimized model, all grid nodes are traversed, and areas with three or more consecutive nodes without elevation values ​​are identified as data gaps. These gaps are filled using the elevation values ​​of valid grid nodes within a 5-meter radius, ensuring that the elevation difference between adjacent nodes does not exceed 0.3 meters to maintain terrain continuity, ultimately generating a complete digital terrain surface model.

[0123] In this embodiment of the invention, by adaptively adjusting the sampling density and interpolation parameters according to the terrain complexity, the range of values ​​is adapted to the terrain diversity of composite rock and soil slopes. This allows the elevation accuracy of the digital terrain surface model to be controlled within 0.1 meters. While preserving key terrain features such as ridges and valleys, interpolation noise is eliminated, providing high-precision terrain benchmark data for subsequent slope displacement vector field calculations and improving the reliability of deformation monitoring under complex terrain conditions.

[0124] In a preferred embodiment of the present invention, based on a digital terrain surface model, iterative optimization is performed by spatial similarity analysis and fusion with a gradient descent algorithm to calculate the slope surface displacement vector field, identify potential sliding surfaces and abnormal deformation areas, and generate displacement field solution results, including:

[0125] The digital terrain surface model is divided into several analysis units. An elevation residual objective function is established within each analysis unit to generate a unit difference quantification dataset. Specifically, when dividing the digital terrain surface model into analysis units, an adaptive gridding method is adopted based on the slope topography complexity. Small analysis units of 5m × 5m are set at the boundary between rock and soil areas and in areas with abrupt slope changes, while large analysis units of 10m × 10m are set in relatively flat and homogeneous areas, ensuring that the unit boundaries are consistent with the direction of the terrain feature lines. For each analysis unit, corresponding areas of the baseline and monitoring periods of the digital terrain surface model are selected, and the elevation differences of all corresponding points within the unit are calculated. The sum of the squares of the elevation differences is taken as the unit elevation residual objective function value to generate a unit difference quantification dataset. The accuracy of the elevation difference calculation is controlled within 0.01 meters.

[0126] Based on the unit differential quantization dataset, an iterative optimization solution is performed using the gradient descent algorithm. This involves calculating the gradient of the objective function and progressively adjusting the spatial position parameters of the analysis units along the negative gradient direction to generate a unit displacement parameter set. Specifically, when using the gradient descent algorithm for iterative optimization based on the unit differential quantization dataset, the initial displacement parameters are set as zero vectors, the learning rate ranges from 0.01 to 0.1, and the maximum number of iterations is set to 50 to 100. During each iteration, the gradient of the objective function of the current unit elevation residual with respect to the displacement parameters is calculated. The gradient direction is the direction in which the function value increases. The displacement parameters are adjusted along the opposite direction (the negative gradient direction), with the adjustment magnitude being the product of the learning rate and the gradient value. After each parameter adjustment, the unit elevation residual objective function value is recalculated. If the difference between the function values ​​of two consecutive iterations is less than 0.05 meters, the iteration stops; otherwise, the adjustment continues until the maximum number of iterations is reached, generating the three-dimensional displacement parameters for each analysis unit, forming a unit displacement parameter set.

[0127] Based on the element displacement parameter set, the displacement vectors of all analysis elements are integrated to construct a global displacement vector field. An optimized displacement vector field is generated through statistical consistency testing, and spatial distribution characteristic analysis is performed to generate potential sliding boundary identification results. Specifically, when constructing the global displacement vector field based on the element displacement parameter set, the displacement parameters of each analysis element are assigned to the element center point. The displacement vector at the element boundary is calculated using a neighborhood interpolation method to ensure continuous change in the displacement vectors of adjacent elements. During the statistical consistency testing, the average difference between the displacement vector of each element and the displacement vectors of 3 to 5 adjacent elements is calculated. If the difference is greater than twice the standard deviation of all element displacement vectors, it is marked as an abnormal element. The displacement vector of the abnormal element is corrected using the weighted average of the surrounding elements, generating an optimized displacement vector field. During the spatial distribution characteristic analysis, the magnitude and direction of the displacement vector of each element are statistically analyzed. Regions with a displacement vector magnitude greater than 0.1 meters or a change in direction angle greater than 30 degrees are marked as potential sliding boundaries, generating potential sliding boundary identification results.

[0128] When implementing the neighborhood interpolation method, the area to be interpolated is divided into grids with a grid side length of 0.5 meters to 2 meters. The grid side length for rock slope areas is 0.5 meters to 1 meter, and the grid side length for soil slope areas is 1 meter to 2 meters. Each grid node is used as an interpolation point. For each interpolation point, a neighborhood is determined centered on that point. The neighborhood radius is adaptively adjusted according to the terrain complexity: 10 to 15 meters for flat terrain, 5 to 10 meters for complex terrain, and 7 to 12 meters for moderate slope terrain. All sample points within the neighborhood are extracted from the structured terrain point set. The horizontal distance between each sample point and the interpolation point is calculated. Valid sample points are selected using a distance threshold, and sample points with a distance greater than the neighborhood radius are removed to ensure that the number of sample points participating in the interpolation is controlled between 10 and 30. Outlier detection is performed on the selected sample points. The mean and standard deviation of the sample point elevation are calculated. Sample points whose elevation values ​​deviate from the mean by more than twice the standard deviation are marked as outliers and removed. The remaining sample points are retained as valid interpolation samples.

[0129] The weighting coefficients of valid sample points are calculated using an inverse distance weighting method. The weight of each sample point is inversely proportional to its horizontal distance to the point to be interpolated; the closer the distance, the greater the weight. The weights of all sample points are normalized so that the sum of the weights is 1, avoiding weight imbalance caused by distance differences. The elevation value of each valid sample point is multiplied by its corresponding weight and then summed to obtain the initial interpolated elevation of the point to be interpolated. The initial interpolation results are verified by randomly selecting 5% to 10% of the points in the structured terrain point set as verification points. The error between the actual elevation and the interpolated elevation of the verification points is calculated, and the average and maximum errors are statistically analyzed. If the maximum error exceeds 0.1 meters or the average error exceeds 0.05 meters, the neighborhood radius is adjusted, and the interpolation process is repeated until the error meets the requirements. The interpolated elevations of all grid nodes are integrated to generate continuous terrain surface elevation data, completing the neighborhood interpolation calculation.

[0130] The results of potential sliding boundary identification are analyzed to determine the spatial morphology and deformation evolution characteristics of the sliding surface, and displacement field calculation results are generated. Specifically, when analyzing the results of potential sliding boundary identification, profile lines are extracted along the direction of the potential sliding boundary in conjunction with slope geological survey data. Based on the changes in the direction and magnitude of displacement vectors on the profile lines, the dip and dip angle of the sliding surface are determined. The dip angle of the sliding surface of rock slopes is taken with reference to the dip of the rock strata, and the dip angle of the sliding surface of soil slopes is taken with reference to the slope angle. When analyzing the deformation evolution characteristics, the changes in the displacement vector field during different monitoring periods are compared, the displacement rate is calculated, and areas with a rate greater than 0.01 m / day are marked as active deformation areas. The deformation is judged by combining the topographic stability index to determine whether the deformation is accelerating. Finally, displacement field calculation results containing the spatial coordinates of the sliding surface, deformation direction, rate, and distribution of active areas are generated.

[0131] In this embodiment of the invention, by using adaptive unit partitioning and gradient descent optimization, the value range adapts to the differences in the terrain of the composite slope, and the displacement calculation accuracy can be controlled within 0.05 meters, improving the accuracy of potential sliding surface identification, providing quantitative basis for slope instability early warning, and effectively ensuring the monitoring reliability of high-risk slopes.

[0132] In a preferred embodiment of the present invention, the displacement field calculation results are input into a risk assessment model, and a quantitative evaluation of slope stability is performed through a multi-source data fusion analysis platform, including:

[0133] Based on the displacement field calculation results, the displacement vector field data and potential sliding surface information contained therein are extracted to generate an initial risk assessment input set. Based on the initial risk assessment input set, the geological survey data and geotechnical parameters of the slope area are integrated to generate a geomechanical parameter set. Specifically, this includes: extracting displacement vector field data and potential sliding surface information from the displacement field calculation results to generate the initial risk assessment input set; for a 500-meter-long slope of a mountain highway from K12+300 to K12+800, a displacement monitoring sampling point grid is determined at a density of 2m×2m, with a total of 12,500 sampling points; based on the displacement field calculation results, three-dimensional displacement coordinate data of each sampling point are obtained for 30 consecutive days at 1-hour intervals, and the values ​​of each sampling point are calculated. The displacement components of the horizontal X-axis and Y-axis and the vertical Z-axis are calculated by dividing the displacement difference by the time interval. The displacement rate ranges from 0.1 to 10 mm / day. The displacement direction angle is calculated using trigonometric functions, with a range of 0 to 360°. These are integrated to form a displacement vector field data. Considering the stratigraphic characteristics of the upper part of the slope, which consists of moderately weathered sandstone interbedded with mudstone, and the lower part, which consists of Quaternary residual slope deposits and silty clay, a displacement mutation threshold is set. The threshold for the sandstone interbedded with mudstone area is 5 mm / day, and the threshold for the silty clay area is 8 mm / day. The displacement gradient threshold is set to 0.02 mm / m. The areas where the displacement exceeds the corresponding area mutation threshold or the displacement gradient exceeds 0.02 mm / m are marked as potential sliding suspicious areas.

[0134] Three-dimensional surface fitting is performed on the displacement data of sampling points in each suspected area to obtain spatial geometric parameters such as the strike, dip angle and dip direction of the potential sliding surface. The strike ranges from 0 to 360°, the dip angle ranges from 30 to 60°, and the dip direction ranges from 0 to 360°. The data is integrated in the format of sampling point number ~ displacement component ~ displacement rate ~ displacement direction angle ~ potential sliding surface number ~ sliding surface geometric parameters to generate the initial risk assessment input set.

[0135] Based on the initial risk assessment input set, geological survey data and geotechnical parameters of the slope area were integrated to generate a geotechnical parameter set. Geological survey data for the slope was collected, including the distribution range of the upper moderately weathered sandstone interbedded with mudstone strata from the top of the slope (K12+300~K12+800) to a height of 60 meters; lithological classification data: sandstone content 60-70%, mudstone content 30-40%; fault and joint development with spacing ranging from 0.5 to 2 meters; stratum thickness: sandstone interbedded with mudstone 50-60 meters; silty clay thickness 5-10 meters. Geotechnical parameters were obtained through indoor direct shear tests and triaxial compression tests. The parameters for sandstone interbedded with mudstone were: cohesion 25-35 kPa; internal friction angle 28-35°; elastic modulus 15-25 GPa; Poisson's ratio 0.2-0.25; and unit weight 22-24 kN / m³. 3The cohesion of silty clay is taken as 15–20 kPa, the internal friction angle as 18–22°, the elastic modulus as 0.5–1.5 GPa, the Poisson's ratio as 0.3–0.35, and the unit weight as 19–21 kN / m³. Based on the strike, dip, and dip of the potential sliding surfaces in the initial input set, the stratigraphic range traversed by each sliding surface is determined. For example, bedding-parallel sliding surfaces mainly traverse sandstone interbedded with mudstone strata, while slump sliding surfaces mainly traverse silty clay strata. The extracted geological survey data of the corresponding strata and geotechnical parameters are associated with the potential sliding surface numbers, resulting in 10 potential sliding surfaces. The Grubbs test with a level of 0.05 is used to detect outliers in the associated data, and the detected outliers are corrected using inverse distance weighted interpolation with a correction radius of 5 meters. The data are organized according to the classification method of potential sliding surfaces to generate a geomechanical parameter set.

[0136] A multi-source fusion database is generated by integrating geomechanical parameter sets with real-time hydrological and meteorological data. Based on this database, the safety factor of each potential sliding surface is calculated, and preliminary stability assessment results are generated. Specifically, this includes: integrating geomechanical parameter sets with real-time hydrological and meteorological data to generate the multi-source fusion database; determining the hydrological and meteorological monitoring indicators for the slope, including groundwater level, pore water pressure, rainfall, atmospheric temperature, and relative humidity. The groundwater level ranges from 100 to 185 meters, with the reference sea level in the area as the zero point; the pore water pressure ranges from 0 to 50 kPa; the rainfall ranges from 0 to 200 mm / day; the atmospheric temperature ranges from -5 to 35℃; and the relative humidity ranges from 30 to 100%. Three sets of groundwater level gauges and pore water pressure gauges are installed at the top, middle, and toe of the slope, and two sets of rain gauges and thermometers / hygrometers are installed at the top. Real-time data are collected at a frequency of once per hour for groundwater level and pore water pressure, once every 5 minutes for rainfall, and once every 30 minutes for temperature and humidity.

[0137] Real-time data is preprocessed, and missing values ​​are filled using linear interpolation, with a single missing value duration not exceeding 2 hours. Noise is filtered using a moving average method with a window size of 5 data points. Based on the spatial coordinates of each potential slip surface in the geomechanical parameter set, a 50-meter radius monitoring area is determined for each slip surface. The geomechanical parameters are spatiotemporally aligned with the hydrological and meteorological data of the corresponding area using second-level timestamps. A data association index is established using the potential slip surface number ~ monitoring time ~ stratigraphic number as index keywords. Data units are unified to international standard units such as pressure unit kPa and length unit m, and the storage format is unified to JSON. The standardized data is stored in a MySQL relational database to generate a multi-source fusion database.

[0138] Safety factors for each potential sliding surface were calculated based on a multi-source fusion database to generate preliminary stability assessment results. The Bishop method was selected as the safety factor calculation method. Based on the strike, dip, and dip direction of each potential sliding surface in the multi-source fusion database, the cohesion, internal friction angle, unit weight, and pore water pressure data of the corresponding strata in the database were substituted into the calculation model. Potential sliding surfaces were divided into several calculation blocks along the strike direction with lengths of 5–8 meters, and each sliding surface was divided into 8–15 blocks. The self-weight of each block was calculated by multiplying the block volume by its unit weight. The seepage force generated by pore water pressure was calculated by multiplying the pore water pressure by the block area. Based on the block geometry and stress state, the normal stress and tangential stress on the sliding surface were calculated, where the normal stress is the normal component of the self-weight minus the pore water pressure, and the tangential stress is the tangential component of the self-weight. The tangential component of the penetrating force is added; based on the Mohr-Coulomb strength criterion, the anti-slip force is calculated by multiplying the cohesion and normal stress by adding the product of the normal stress and the tangent of the internal friction angle, and the sliding force is calculated by the tangential stress. The safety factor of the entire sliding surface is solved by an iterative method, with the number of iterations not exceeding 20 and the convergence error less than 0.001. A safety factor evaluation threshold is set, where a safety factor ≥ 1.25 indicates a stable state, 1.05 ≤ safety factor < 1.25 indicates a basically stable state, 1.0 ≤ safety factor < 1.05 indicates a sub-stable state, and a safety factor < 1.0 indicates an unstable state. Based on the comparison between the calculated safety factor value of each sliding surface and the threshold, the calculated safety factor value ranges from 0.8 to 1.5 to determine the stability state. The data is organized according to the format of the calculation parameters for the potential sliding surface number to generate preliminary stability assessment results.

[0139] Based on the preliminary stability assessment results, stress-strain analysis results are generated. The stress-strain analysis results are then fused with the preliminary stability assessment results, and a comprehensive stability index is calculated using a weighted fusion algorithm to generate a slope stability level evaluation set. Specifically, this includes: generating stress-strain analysis results based on the preliminary stability assessment results; conducting stress-strain analysis using a combination of on-site monitoring and theoretical calculations; and deploying stress sensors at potential sliding surface areas and stratigraphic boundaries in the slope based on the lithological data and potential sliding surface information from a multi-source fusion database. The sensor deployment density is [missing information]. One sensor is installed per 100 square meters, with a total of 50 sensors deployed. The sensor monitoring area is determined using tetrahedral element meshing logic. The monitoring area size for potential sliding surfaces and stratigraphic boundaries is set to 3m×3m×3m, while other areas are set to 5m×5m×5m. The elastic modulus and Poisson's density of each stratum in the database are used as the theoretical calculation basis parameters. For sliding surface areas initially assessed as underlying or unstable, the frictional characteristic parameters of the monitoring area are determined according to the corresponding internal friction angle tangent value. The frictional characteristic parameter values ​​range from 0.33 to 0.70. Simultaneously, the reference value for normal force is determined to be 5×10⁻⁶. 5 ~1×10 6 N / m³.

[0140] During monitoring, the boundary stress state at each sensor location was recorded. Displacement in the X, Y, and Z axes was restricted at the bottom, and in the X and Y axes at the sides. There was no displacement restriction at the top. The uniformly distributed load magnitude for each monitoring area was determined based on real-time pore water pressure data from the database. Stress data of the slope area was collected in real-time by sensors, and strain data was derived through theoretical calculations to obtain the distribution of maximum principal stress, minimum principal stress, shear stress, axial strain, shear strain, and volumetric strain. The maximum principal stress ranged from 0 to 200 kPa, the minimum principal stress from ~50 to 100 kPa, the shear stress from 0 to 80 kPa, and the axial strain from ~5 × 10⁻⁵ kPa. -4 ~5×10 -4 The shear strain ranges from 0 to 3 × 10⁻⁶. -4 The volumetric strain range is ~3×10 -4 ~3×10 -4 Extract the maximum shear stress value, stress concentration factor, and maximum shear strain value of each potential sliding surface region. The stress concentration factor is the ratio of the maximum shear stress to the average shear stress, with a value range of 1.2 to 2.5. Plot the stress-strain distribution cloud map, integrate the characteristic parameters and the distribution cloud map, and generate the stress-strain analysis results.

[0141] The results of stress-strain analysis and preliminary stability assessment were combined using a weighted fusion algorithm to calculate a comprehensive stability index, generating a slope stability level evaluation set. Four evaluation indicators were determined: the safety factor A from the preliminary assessment results, the maximum shear stress value B from the stress-strain analysis results, the stress concentration factor C, and the maximum shear strain value D. The weights were determined using the analytic hierarchy process (AHP), a judgment matrix was constructed, and eigenvectors were calculated using the sum-product method. After a consistency check, the weights were determined, with a consistency ratio CR < 0.1. The weights were A = 0.5, B = 0.2, C = 0.2, and D = 0.1. The indicators are standardized. The safety factor is positively standardized, calculated as the difference between the actual value and 1.0, divided by the difference between 1.5 and 1.0, with a range of 0 to 1. The maximum shear stress is negatively standardized, calculated as the difference between 80 and the actual value, divided by the difference between 80 and 0, with a range of 0 to 1. The stress concentration factor is negatively standardized, calculated as the difference between 2.5 and the actual value, divided by the difference between 2.5 and 1.2, with a range of 0 to 1. The maximum shear strain is negatively standardized, calculated as 3 × 10⁻⁶. -4 Difference between the actual value and the actual value, divided by 3 × 10 -4Subtract 0 from the difference, with a value ranging from 0 to 1; calculate the index value of each potential sliding surface according to the formula A×0.5+B×0.2+C×0.2+D×0.1, with the comprehensive stability index value ranging from 0 to 1; set the grade range: comprehensive stability index ≥ 0.8 indicates Level 1 stability; 0.6≤comprehensive stability index<0.8 indicates Level 2 basic stability; 0.4≤comprehensive stability index<0.6 indicates Level 3 understability; comprehensive stability index<0.4 indicates Level 4 instability; determine the stability grade of each sliding surface by comparing the calculated value with the grade range, organize the data according to the format of potential sliding surface number, and generate a slope stability grade evaluation set.

[0142] Based on the slope stability level evaluation set, slope risk areas are delineated and early warning levels are determined, generating slope risk zoning and early warning schemes. Based on the slope risk zoning and early warning schemes, a quantitative slope stability evaluation report is generated, specifically including: delineating slope risk areas and determining early warning levels based on the slope stability level evaluation set to generate slope risk zoning and early warning schemes; dividing the actual length of the slope from K12+300 to K12+800 into 10 evaluation units of 50 meters each, with each unit corresponding to K12+300 to K12+350 to K12+750 to K12+800, and each unit containing 1 to 2 potential sliding surfaces; determining the dominant stability level of each evaluation unit: when a level 4 instability sliding surface exists within the unit, the dominant level is level 4; if there is no level 4 but a level 3 surface exists, the dominant level is level 3, and so on; dividing risk areas according to the dominant level: level 1 stability corresponds to low-risk areas, level 2 basically stable corresponds to relatively low-risk areas, level 3 understability corresponds to medium-risk areas, and level 4 instability corresponds to high-risk areas.

[0143] Based on displacement rate data from the displacement vector field in the multi-source fusion database, early warning level standards are set: blue warning for low-risk areas with a displacement rate <0.5 mm / day; yellow warning for lower-risk areas or 0.5 mm / day ≤ displacement rate <1.0 mm / day; orange warning for medium-risk areas or 1.0 mm / day ≤ displacement rate <2.0 mm / day; and red warning for high-risk areas or displacement rate ≥2.0 mm / day. The warning level is determined according to the dominant level and displacement rate of each evaluation unit, with the displacement rate ranging from 0.1 to 5.0 mm / day. A spatial distribution map of the risk area is drawn using GIS software, marking the range of each area, such as K12+400~K12+450 as a medium-risk area, dominant level, warning level, and judgment criteria, such as a safety factor of 1.03 and a displacement rate of 1.2 mm / day. This integration forms a slope risk zoning and early warning scheme. The beneficial effects of this process are as follows: dividing the evaluation units according to the actual road sections ensures that the zoning is consistent with the actual project; the dual judgment of dominant level and displacement rate improves the accuracy of early warning; GIS visualization facilitates intuitive understanding of risk distribution; it solves the problem of vague early warning in traditional risk zoning; and it provides a clear target for disaster prevention and control of high-risk slopes with houses on the top and highways at the foot.

[0144] A quantitative evaluation report on slope stability was generated based on the slope risk zoning and early warning scheme. The report's structural framework was determined, including six parts: project overview, data sources and processing, detailed evaluation steps, evaluation results, risk prevention and control recommendations, and conclusions. In the project overview section, the geographical location of the slope was described in detail as the K12+300~K12+800 section of a mountainous expressway; its scale was 500 meters long, with a maximum slope height of 85 meters and an average slope of 55°; its stratigraphic characteristics were moderately weathered sandstone interbedded with mudstone in the upper part and Quaternary residual colluvial silty clay in the lower part; its surrounding environment included residential houses at the top of the slope and the main line of the expressway at the foot of the slope; and its risk characteristics were bedding-parallel sliding and potential landslides during the rainy season. In the data sources and processing section, it was explained that the displacement field calculation data came from the slope monitoring system using a 2m×2m sampling grid; the geological survey data came from borehole surveys with a total of 8 boreholes; and the hydrological and meteorological data came from on-site monitoring equipment, including 3 sets of groundwater level gauges, etc. The preprocessing methods for each data were also explained, such as linear interpolation and moving average.

[0145] The detailed process of each evaluation step is described sequentially, outlining the generation process of the initial input set, geomechanical parameter set, multi-source fusion database, preliminary assessment results, stress-strain analysis results, grade evaluation set, and risk zoning scheme. The methods used in each step are clearly defined, such as Kriging interpolation, Bishop's method, and analytic hierarchy process; parameter values, such as cohesion of 25–35 kPa, safety factor threshold of 1.05; and calculation details, such as block division of 5–8 meters. The evaluation results section presents the safety factors for each potential sliding surface, such as 1.12 for sliding surface 1 and 0.98 for sliding surface 5, and the comprehensive stability index, such as 0.75 for sliding surface 1 and 0.75 for sliding surface 5. 0.38. Stability level and risk area distribution, such as the K12+400~K12+500 section being a high-risk area, and the warning level, such as the red warning for this section, and the corresponding GIS distribution map; In the risk prevention and control recommendations section, for the red warning area, such as the K12+400~K12+500 section, it is proposed to set up anchor support anti-slide piles and optimize the drainage system. Among them, the anchor length is 8~10 meters with a spacing of 2m×2m, the pile length is 15~20 meters, the pile diameter is 1.2 meters, and the spacing is 5 meters. The drainage system optimization includes adding intercepting ditches and blind ditches. For the orange warning area, it is recommended to increase the monitoring frequency to 15 minutes / time.

[0146] For yellow and blue alert areas, recommendations were made for routine monitoring at 1 hour / time and regular inspection at 1 week, respectively. The conclusion summarized the overall stability of the slope, e.g., 30% of the area was high-risk and 20% was medium-risk. Key risk areas, such as the K12+400~K12+500 section, and priority control areas were identified. The evaluation results were explained as applicable under the condition that there were no significant changes in the strata or extreme meteorological anomalies in the past 3 months, and were limited by the fact that the impact of seismic loads was not considered. The report data was reviewed and compared three times, with an error of less than 5%. Logical verification ensured consistency in data transmission at each step, and the terminology conformed to geotechnical engineering investigation standards, resulting in a complete quantitative evaluation report of slope stability.

[0147] In this embodiment of the invention, the segmentation and iterative calculation ensure the accuracy of the safety factor calculation. Combined with real-time pore water pressure data, it can dynamically reflect the impact of the rainy season on stability. Clear threshold division can quickly determine the risk status of the sliding surface, solving the problem of low efficiency and inability to adapt to dynamic environments in traditional methods. Standardization eliminates differences in index dimensions, the weights determined by the analytic hierarchy process conform to the actual focus of the project, and clear level division facilitates intuitive judgment of the degree of risk, improving the comprehensiveness and reliability of the evaluation.

[0148] Embodiments of the present invention also provide a computing device, including: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the system as described above. All implementations in the above system embodiments are applicable to this embodiment and can achieve the same technical effects.

[0149] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the system as described above. All implementations in the above system embodiments are applicable to this embodiment and can achieve the same technical effects.

[0150] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A slope displacement monitoring data processing system based on unmanned aerial vehicle laser radar, characterized in that, The method comprises the following steps: a collection module is used to acquire an original laser radar observation data set; a processing module is used to preprocess the original laser radar observation data set to obtain an optimized laser radar observation data set; a calculation module is used to extract stable terrain feature points in the observation data, filter out the ground protrusions, rock layer junction lines and terrain abrupt change points as high-reliability feature points based on the curvature and normal vector change rate indexes, generate an initial feature point set, construct a feature descriptor for the initial feature point set, generate a feature description vector based on the local surface geometric properties, obtain a standardized feature description set, perform similarity matching on the standardized feature description set, eliminate the mis-matching points through bidirectional consistency inspection, generate a reliable feature point corresponding relationship set, calculate the final spatial transformation parameters based on the reliable feature point corresponding relationship set, and generate a coordinate transformation parameter set; the coordinate transformation parameter set is used to perform rigid transformation on all the laser radar observation data, convert the laser radar observation data into a reference coordinate system, and generate a spatial reference unified data set; the spatial reference unified data set is fused to eliminate the redundant points in the data overlapping area, and a complete slope three-dimensional observation data set is generated; a generation module is used to generate a directional reference framework based on the slope three-dimensional observation data set, and generate an adaptive correction factor according to the topological relationship and structural features of each region in the framework, comprising the following steps: a first spatial reference point located in the stable bedrock region at the top of the slope and a second spatial reference point located in the stable bedrock region at the toe of the slope are selected from the slope three-dimensional observation data set to generate an initial reference axis; the initial reference axis is used as a reference direction to construct a regional coordinate system with the first spatial reference point as the origin and the reference line direction as the main axis, and a directional reference framework is generated; based on the directional reference framework, the monitoring area is divided into a plurality of independent analysis units according to the slope geological structure features and the terrain change trend, and a unit division result is generated; the terrain structure properties of each analysis unit are extracted to generate a unit terrain feature description set; based on the relative spatial position relationship of each analysis unit and the two spatial reference points, the azimuth weight and distance weight of each unit relative to the reference line are calculated to generate a unit weight distribution set; the unit terrain feature description set and the unit weight distribution set are fused to calculate the terrain stability comprehensive index of each analysis unit to generate a unit stability evaluation set; based on the unit stability evaluation set, the stability index is converted into a ground object classification parameter adjustment coefficient to generate an adaptive correction factor; a classification module optimizes the ground object classification processing parameters according to the adaptive correction factor to realize accurate separation of vegetation coverage, artificial structures and natural ground surfaces, and generates an optimized terrain feature data set; a construction module is used to perform spatial interpolation processing on the terrain feature data set to construct a digital terrain surface model; an analysis module is used to calculate a slope surface displacement vector field based on the digital terrain surface model, identify potential sliding surfaces and deformation abnormal areas, and generate a displacement field solution result through spatial similarity analysis and iterative optimization by fusing the gradient descent algorithm. The output module is configured to input the displacement field calculation result into a risk assessment model, and perform quantitative evaluation on the slope stability through a multi-source data fusion analysis platform. 2.The UAV LiDAR-based slope displacement monitoring data processing system of claim 1, wherein, The original laser radar observation data set is preprocessed to obtain an optimized laser radar observation data set, including: The original laser radar observation data set is subjected to noise filtering, and discrete noise points are removed based on statistical outlier analysis to generate a preliminary purified data set; the preliminary purified data set is subjected to intensity correction, and the intensity distortion caused by the measurement geometric effect is eliminated through distance and incident angle normalization processing to generate an intensity normalized data set; The intensity normalized data set is subjected to coordinate correction, and the geometric deviation caused by the attitude change of the unmanned aerial vehicle platform is compensated based on the sensor pose parameters and inertial measurement unit data to generate a geometric corrected data set; The geometric corrected data set is subjected to sampling optimization, and the point density is adaptively adjusted based on the slope terrain features to remove redundant data points to generate the optimized laser radar observation data set. 3.The UAV LiDAR-based slope displacement monitoring data processing system of claim 2, wherein, The ground feature classification processing parameters are optimized according to the adaptive correction factor to realize accurate separation of vegetation coverage, artificial structures and natural ground surfaces, and an optimized terrain feature data set is generated, including: According to the adaptive correction factor, the analysis unit parameter adjustment coefficient is analyzed to generate a unit parameter adjustment set; based on the unit parameter adjustment set, the classification threshold parameter in the ground feature classification algorithm is adaptively adjusted according to the terrain stability features of different analysis units to generate an optimized classification parameter set; According to the optimized classification parameter set, the initial ground feature classification result is generated; the initial classification result is processed to generate an optimized classification result; Based on the optimized classification result, the spatial distribution information of the vegetation coverage area and the artificial structure area is extracted to generate a ground feature mask data set; the ground feature mask data set is used to process the slope stereoscopic observation data set to remove the data points corresponding to the vegetation and artificial structures to generate a bare ground three-dimensional point set; The bare ground three-dimensional point set is processed to generate an optimized terrain feature data set. 4.The UAV LiDAR-based slope displacement monitoring data processing system of claim 3, wherein, The terrain feature data set is subjected to spatial interpolation processing to construct a digital terrain surface model, including: The terrain feature data set is subjected to terrain feature analysis, and the terrain feature points are identified based on the local curvature change and the elevation variation coefficient; the ridge lines, valley lines and slope mutation points are screened as key terrain points to generate a candidate terrain feature point set; the candidate terrain feature point set is subjected to density optimization sampling, and the sampling density is adaptively adjusted based on the terrain complexity to sparsely sample in flat terrain areas and densely sample in complex areas to generate a structured terrain point set; Based on the structured terrain point set, a spatial interpolation function is constructed to obtain final interpolation parameters to generate an optimized interpolation function; the optimized interpolation function is used to calculate the numerical values of regularly distributed elevation points to generate a preliminary terrain surface model; The preliminary terrain surface model is subjected to smoothing processing to eliminate local interpolation anomalies to generate an optimized terrain surface model; the optimized terrain surface model is subjected to integrity inspection, and the data vacancy areas are filled based on the terrain continuity principle to generate a final digital terrain surface model. 5.The UAV LiDAR-based slope displacement monitoring data processing system of claim 4, wherein, Based on the digital terrain surface model, the spatial similarity analysis is performed and the gradient descent algorithm is fused for iterative optimization to calculate the slope surface displacement vector field, identify the potential sliding surface and deformation anomaly area, and generate the displacement field calculation result, including: The digital terrain surface model is divided into a plurality of analysis units, and an elevation residual objective function is established in each analysis unit to generate a unit difference quantization data set; Based on the unit difference quantization data set, the gradient descent algorithm is used for iterative optimization to calculate the gradient of the objective function and gradually adjust the spatial position parameters of the analysis unit along the negative gradient direction to generate a unit displacement parameter set; Based on the unit displacement parameter set, the displacement vectors of all analysis units are integrated to construct a global displacement vector field, and the optimization displacement vector field is generated through statistical consistency test processing, and the displacement spatial distribution characteristic analysis is performed to generate the potential sliding boundary identification result; The potential sliding boundary identification result is analyzed to determine the spatial form and deformation evolution characteristics of the sliding surface, and the displacement field calculation result is generated. 6.The UAV LiDAR-based slope displacement monitoring data processing system of claim 5, wherein, The displacement field calculation result is input into the risk assessment model, and the slope stability is quantitatively evaluated through the multi-source data fusion analysis platform, including: According to the displacement field calculation result, the displacement vector field data and potential sliding surface information contained therein are extracted to generate an initial risk assessment input set; based on the initial risk assessment input set, the geological survey data and geotechnical mechanics parameters of the slope area are integrated to generate a geomechanics parameter set; Fusion of geomechanics parameter set and real-time monitoring of hydrological and meteorological data, generate multi-source fusion database; based on the multi-source fusion database, calculate the safety factor of each potential sliding surface to generate the preliminary stability evaluation result; Based on the preliminary stability evaluation result, the stress-strain analysis result is generated; the stress-strain analysis result and the preliminary stability evaluation result are fused, and the comprehensive stability index is calculated by weighted fusion algorithm to generate the slope stability grade evaluation set; Based on the slope stability grade evaluation set, the slope risk area is divided and the warning level is determined to generate the slope risk zoning and warning scheme; according to the slope risk zoning and warning scheme, the slope stability quantitative evaluation report is generated. 7.The UAV LiDAR-based slope displacement monitoring data processing system of claim 6, wherein, Based on the unit difference quantization data set, the gradient descent algorithm is used for iterative optimization to calculate the gradient of the objective function and gradually adjust the spatial position parameters of the analysis unit along the negative gradient direction to generate a unit displacement parameter set, including: Initialize the displacement parameters of each analysis unit to generate an initial displacement parameter set, and calculate the elevation residual objective function value of each analysis unit to generate an initial objective function value set; Iterative optimization is performed on the initial displacement parameter set, and in each iteration, the partial derivative of the objective function with respect to the displacement parameter of each analysis unit is calculated to generate a gradient vector set, and the update direction and step size of each analysis unit displacement parameter are determined, and the displacement parameter is adjusted along the negative gradient direction to generate an updated displacement parameter set; The updated displacement parameter set is recalculated to generate a new objective function value set; compare the new objective function value set with the objective function value set of the last iteration to determine whether the iteration meets the preset convergence condition, and finally generate the unit displacement parameter set.

8. A computing device, comprising: including: one or more processors; a storage device for storing one or more programs, when executed by the one or more processors, cause the one or more processors to implement a system as claimed in any of claims 1 to 7.

Images