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

Abstract

The invention provides a slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar, and relates to the technical field of data processing, and the system comprises the steps: carrying out the spatial interpolation processing of a topographic feature data set, and constructing a digital topographic surface model; the displacement field calculation module is used for performing iterative optimization based on a digital terrain surface model through spatial similarity analysis and fusion with a gradient descent algorithm, calculating a slope surface displacement vector field, identifying a potential sliding surface and a deformation abnormal region, and generating a displacement field calculation result; and the evaluation module is used for inputting a displacement field calculation result into a risk evaluation model and carrying out slope stability quantitative evaluation through a multi-source data fusion analysis platform. According to the invention, the practicability and operability of the monitoring result are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of data processing, in particular to a slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar. BACKGROUND

[0002] Slope stability monitoring is a core technical problem to ensure the safety of engineering in the specific section of mountainous highway slope.

[0003] The existing data processing procedure has not formed a standardized system, and the monitoring accuracy and efficiency in complex environment cannot meet the engineering requirements. For example, the existing monitoring technology relies on inclinometer monitoring or osmometer monitoring, and such methods have the following technical limitations: the slope is long and high, the inclinometer and osmometer are single-point arranged (usually a certain distance apart to arrange monitoring sections, and a small number of monitoring points are arranged at each section), and only the local point deep deformation and osmotic pressure can be reflected. For example, within a certain thickness range of lower silty clay, if the local osmotic pressure of a certain section of soil increases suddenly due to uneven distribution of block stones, a single-point osmometer cannot capture the anomaly, and early warning of collapse is easy to miss. The lower silty clay is loose and contains block stones, and the inclinometer needs to be drilled to the deep layer. During the drilling process, the original structure of the soil is easily damaged, which reduces the local soil shear strength and may even induce small-scale collapse, which is contrary to the original intention of monitoring and protecting the slope. SUMMARY

[0004] The technical problem to be solved by the present application is to provide a slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar, which improves the practicality and operability of the monitoring results.

[0005] To solve the above technical problems, the technical scheme of the present application is as follows:

[0006] In a first aspect, the slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar comprises:

[0007] The slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar comprises:

[0008] The slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar comprises:

[0009] The acquisition module is used to obtain an original laser radar observation data set;

[0010] The processing module is used to preprocess the original laser radar observation data set to obtain an optimized laser radar observation data set;

[0011] The computing module is configured to extract stable terrain feature points in observation data, screen out ground surface protrusions, rock layer junction lines and terrain abrupt change points as high-reliability feature points based on a curvature and normal vector change rate index, and generate an initial feature point set; perform feature descriptor construction on the initial feature point set, generate a feature description vector based on local surface geometric properties, and obtain a standardized feature description set; perform similarity matching on the standardized feature description set, eliminate mis-matching point pairs through bidirectional consistency inspection, and generate a reliable feature point corresponding relationship set; calculate final spatial transformation parameters based on the reliable feature point corresponding relationship set, and generate a coordinate transformation parameter set; perform rigid transformation on all laser radar observation data using the coordinate transformation parameter set, uniformly convert the laser radar observation data to a reference coordinate system, and generate a spatial reference unified data set; and perform fusion on the spatial reference unified data set, eliminate redundant points in a data overlapping area, and generate a complete slope three-dimensional observation data set.

[0012] The generating module is configured to establish a directional reference framework based on the slope three-dimensional observation data set, and generate an adaptive correction factor according to a topological relationship and a structural feature of each area in the framework.

[0013] The classification module is configured to optimize ground object classification processing parameters according to the adaptive correction factor, realize accurate separation of vegetation coverage, artificial structures and natural ground surfaces, and generate an optimized terrain feature data set.

[0014] The constructing module is configured to perform spatial interpolation processing on the terrain feature data set, and construct a digital terrain surface model.

[0015] The analyzing module is configured to calculate a slope surface displacement vector field, identify a potential sliding surface and a deformation abnormal area, and generate a displacement field solution result based on the digital terrain surface model, through spatial similarity analysis and iterative optimization of a gradient descent algorithm.

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

[0017] Further, the original laser radar observation data set is preprocessed to obtain an optimized laser radar observation data set, including:

[0018] The original laser radar observation data set is subjected to noise filtering, discrete noise points are removed based on statistical outlier analysis, and a preliminary purified data set is generated; the preliminary purified data set is subjected to intensity correction, distance and incident angle normalization processing is performed, intensity distortion caused by measurement geometric effects is eliminated, and an intensity normalized data set is generated.

[0019] The intensity normalized data set is coordinate corrected, geometric deviation caused by UAV platform attitude change is compensated based on sensor pose parameters and inertial measurement unit data, and a geometric correction data set is generated;

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

[0021] Further, based on the three-dimensional observation data set of the slope, two stable reference bases are preset, a directional reference frame is established, and adaptive correction factors are generated according to the topological relationship and structural features of each region in the frame, including:

[0022] The first spatial reference point located in the stable bedrock region at the top of the slope and the second spatial reference point located in the stable bedrock region at the toe of the slope are selected from the three-dimensional observation data set of the slope, and an initial reference axis is generated;

[0023] Taking the initial reference axis as the reference direction, a regional coordinate system with the first spatial reference point as the origin and the reference line direction as the main axis is constructed, and a directional reference frame is generated;

[0024] Based on the directional reference frame, the monitoring area is divided into several independent analysis units according to the geological structure features and terrain change trend of the slope, and a unit division result is generated;

[0025] The topographic structure attributes of each analysis unit are extracted to generate a unit topographic 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;

[0026] The unit topographic feature description set and the unit weight distribution set are fused to calculate the topographic stability comprehensive index of each analysis unit to generate a unit stability evaluation set;

[0027] 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.

[0028] Further, the ground object classification processing parameters are optimized according to the adaptive correction factor to realize accurate separation of vegetation cover, artificial structures and natural ground surface, and an optimized topographic feature data set is generated, including:

[0029] According to the adaptive correction factor, the analysis unit parameter adjustment set is generated by analyzing the parameter adjustment coefficient of each analysis unit; based on the unit parameter adjustment set, the classification threshold parameter in the ground object classification algorithm is adaptively adjusted according to the topographic stability features of different analysis units to generate an optimized classification parameter set;

[0030] The preliminary feature classification result is generated by performing preliminary feature classification on the slope stereoscopic observation data set according to the optimized classification parameter set; and the optimized classification result is generated by processing the preliminary feature classification result;

[0031] 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 feature mask data set; and the feature mask data set is used to process the slope stereoscopic observation data set to remove the data points corresponding to the vegetation and the artificial structure, thereby generating a bare ground three-dimensional point set;

[0032] The bare ground three-dimensional point set is processed to generate an optimized terrain feature data set.

[0033] Further, the terrain feature data set is subjected to spatial interpolation processing to construct a digital terrain surface model, including:

[0034] The terrain feature analysis is performed on the optimized terrain feature data set, the terrain feature points are identified based on the local curvature change and the elevation variation coefficient, the ridge lines, the valley lines and the slope mutation points are screened as the key terrain points, and a candidate terrain feature point set is generated; the density optimization sampling is performed on the candidate terrain feature point set, the sampling density is adaptively adjusted based on the terrain complexity, the sampling is sparse in the flat terrain area and dense in the complex area, and a structured terrain point set is generated;

[0035] Based on the structured terrain point set, a spatial interpolation function is constructed to obtain final interpolation parameters, and an optimized interpolation function is generated; the optimized interpolation function is used to calculate the numerical values of the regularly distributed elevation points, and a preliminary terrain surface model is generated;

[0036] The preliminary terrain surface model is subjected to smoothing processing to eliminate local interpolation abnormalities, and an optimized terrain surface model is generated; the optimized terrain surface model is subjected to integrity inspection, and the data vacancy area is filled based on the terrain continuity principle, and a final digital terrain surface model is generated.

[0037] Further, based on the digital terrain surface model, the spatial similarity analysis is performed and the gradient descent algorithm is fused for iterative optimization, the slope surface displacement vector field is calculated, the potential sliding surface and the deformation abnormal area are identified, and a displacement field solution result is generated, including:

[0038] 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;

[0039] Based on the unit difference quantization data set, the gradient descent algorithm is used for iterative optimization solution, the gradient of the objective function is calculated, and the spatial position parameters of the analysis unit are gradually adjusted along the negative gradient direction, thereby generating a unit displacement parameter set;

[0040] Based on the unit displacement parameter set, the displacement vector of all analysis units is integrated to construct a global displacement vector field, and an optimized displacement vector field is generated through statistical consistency test processing, and the spatial distribution characteristics of displacement are analyzed to generate the potential sliding boundary identification result;

[0041] 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.

[0042] Further, the displacement field calculation result is input into the risk assessment model, and the slope stability quantitative evaluation is carried out through the multi-source data fusion analysis platform, including:

[0043] 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;

[0044] The geomechanics parameter set and real-time monitored hydro-meteorological data are fused to generate a multi-source fusion database; based on the multi-source fusion database, the safety factor of each potential sliding surface is calculated to generate a preliminary stability evaluation result;

[0045] Based on the preliminary stability evaluation result, stress-strain analysis results are generated; the stress-strain analysis results and the preliminary stability evaluation results are fused, and a comprehensive stability index is calculated through a weighted fusion algorithm to generate a slope stability grade evaluation set;

[0046] Based on the slope stability grade evaluation set, the slope risk area is divided and the warning level is determined to generate a slope risk zoning and warning scheme; according to the slope risk zoning and warning scheme, a slope stability quantitative evaluation report is generated.

[0047] Further, based on the unit difference quantization data set, a gradient descent algorithm is used for iterative optimization solution, the gradient of the objective function is calculated, and the spatial position parameters of the analysis unit are gradually adjusted along the negative gradient direction to generate the unit displacement parameter set, including:

[0048] The displacement parameters of each analysis unit are initialized to generate an initial displacement parameter set, and the elevation residual objective function value of each analysis unit is calculated to generate an initial objective function value set;

[0049] The initial displacement parameter set is iteratively optimized, and in each iteration, the partial derivative of the objective function of each analysis unit with respect to the displacement parameter is calculated to generate a gradient vector set, and the update direction and step size of the displacement parameter of each analysis unit are determined, and the displacement parameter is adjusted along the negative gradient direction to generate an 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 last iteration to determine whether the iteration meets the preset convergence condition, and finally the unit displacement parameter set is generated.

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

[0052] one or more processors;

[0053] a storage device storing one or more programs, when the one or more programs are executed by the one or more processors, the one or more processors implement the system.

[0054] In a third aspect, a computer-readable storage medium stores a program, when the program is executed by a processor, the system is implemented.

[0055] The above scheme of the present application at least includes the following beneficial effects:

[0056] By adopting the feature point matching algorithm and the space transformation parameter optimization technology, the registration error in the sparse feature area such as smooth bare rock and broken terrain is effectively controlled and reduced compared with the prior art. Through unified space reference conversion and multi-period data fusion, the consistency of cross-time sequence monitoring data is ensured, and reliable data foundation is provided for displacement analysis.

[0057] Through the fine classification algorithm of ground objects, the accurate separation of vegetation coverage, artificial structures and natural ground surface is realized, the classification accuracy in the low-vegetation coverage area is improved, and the distortion problem of terrain feature extraction in the low shrub area is effectively solved. The optimized terrain feature data set provides a high-quality data source for the construction of the digital terrain surface model, and improves the model accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0058] Figure 1 is a schematic diagram of a slope displacement monitoring data processing system based on an unmanned aerial vehicle laser radar provided by an embodiment of the present application.

[0059] Figure 2 is a flowchart of extracting and matching stable terrain feature points in the observation data set, calculating the space transformation parameters by using the feature point matching algorithm, accurately registering the optimized laser radar observation data, completing the conversion and fusion of all data to the unified space reference, and generating the slope stereoscopic observation data set. DETAILED DESCRIPTION

[0060] Exemplary embodiments of the present disclosure will be described in greater detail below with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it is understood that the present disclosure can be embodied in various forms and should not be limited by 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 present disclosure to those skilled in the art.

[0061] As Figure 1 shown, the embodiment of the present application proposes a slope displacement monitoring data processing system based on unmanned aerial vehicle laser radar, comprising:

[0062] The slope displacement monitoring data processing system based on unmanned aerial vehicle laser radar comprises:

[0063] The acquisition module is configured to obtain an original laser radar observation data set.

[0064] The processing module is configured to preprocess the original laser radar observation data set to obtain an optimized laser radar observation data set.

[0065] The computing module is configured to extract and match stable terrain feature points in the observation data set, calculate spatial transformation parameters using a feature point matching algorithm, accurately register the optimized laser radar observation data, complete the conversion and fusion of all data to a unified spatial reference, and generate a slope stereo observation data set.

[0066] The generation module is configured to generate two stable reference bases based on the slope stereo observation data set, establish a directional reference frame, and generate an adaptive correction factor according to the topological relationship and structural features of each region in the frame.

[0067] The classification module optimizes the ground feature classification processing parameters according to the adaptive correction factor, realizes accurate separation of vegetation cover, artificial structures and natural ground surface, and generates an optimized terrain feature data set.

[0068] The construction module is configured to perform spatial interpolation processing on the terrain feature data set to construct a digital terrain surface model.

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

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

[0071] In the embodiment of the present application, by adopting the feature point matching algorithm and the space transformation parameter optimization technology, the registration error in the sparse feature area such as smooth bare rock and broken terrain is effectively controlled, which is lower than the prior art, and through the unified space reference conversion and multi-period data fusion, the consistency of the cross-time sequence monitoring data is ensured, thereby providing a reliable data basis for displacement analysis.

[0072] The accurate separation of vegetation coverage, artificial structures and natural ground is realized through the fine classification algorithm of ground objects, the classification accuracy in the low and medium vegetation coverage area is improved, the distortion problem of terrain feature extraction in the low shrub area is effectively solved, the optimized terrain feature data set provides a high-quality data source for the construction of the digital terrain surface model, and the model accuracy is improved.

[0073] In a preferred embodiment of the present application, the original laser radar observation data set is preprocessed to obtain an optimized laser radar observation data set, including:

[0074] 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, specifically including: when the original laser radar observation data set is subjected to noise filtering, for each three-dimensional point in the data set, all neighborhood points with a preset radius range of 0.3 meters to 1 meter centered on the three-dimensional point are selected, the Euclidean distance between the three-dimensional point and each neighborhood point is calculated, the average distance and the standard deviation are calculated based on all Euclidean distances, the neighborhood points with the Euclidean distance greater than the sum of the average distance and three times the standard deviation are marked as potential noise points, and after removing all potential noise points, the preliminary purified data set is generated. When the preliminary purified data set is subjected to intensity correction, the actual measurement distance between the emission end and each three-dimensional point is obtained through the ranging module of the laser radar sensor, the distance range is 50 meters to 200 meters, the incident angle of the laser beam when reaching the three-dimensional point is obtained through the angle measurement unit built in the sensor, and the incident angle range is 15 degrees to 60 degrees.

[0075] In constructing the intensity normalization model, the correlation between the original intensity value and the actual measured distance and the incident angle is determined, the original intensity value is taken as the dependent variable, the reciprocal of the square of the actual measured distance and the cosine value of the incident angle are taken as the independent variables, and a multiple regression model framework is established; in the model training stage, a standard reference target with a known reflectivity in the slope area is selected, distance measuring points are set at intervals of 10 meters within a distance range of 50 meters to 200 meters, the laser emission angle is adjusted at intervals of 5 degrees within an incident angle range of 15 degrees to 60 degrees at each distance measuring point, the original intensity data of the standard target are collected, and a training sample set is formed; the training sample set is subjected to data fitting, the distance correction coefficient and the angle correction coefficient are calculated, the distance correction coefficient has a value range of 0.8 to 1.2, the angle correction coefficient has a value range of 0.9 to 1.1, the deviation between the model output value and the theoretical reflection intensity of the standard target is controlled to be within 5% through least square method optimization of the coefficients; in the least square method optimization coefficient calculation process, an error function is established, the square of the difference between the model output correction intensity value of each training sample and the theoretical reflection intensity of the standard target is taken 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 added to obtain a total error function, the total error function expression is the sum of each sample error term, and the sample error term is equal to the square of the model output correction intensity value of the standard target theoretical reflection intensity.

[0076] The partial derivatives of the total error function with respect to the distance correction coefficient and the angle correction coefficient are calculated, a binary linear equation set is established according to the extreme condition that the partial derivatives are zero, the equation set contains two unknowns, the distance correction coefficient and the angle correction coefficient, and the coefficients of the equation are calculated from the actual measured distance, the incident angle and the original intensity value in the training sample set; the binary linear equation set is solved by a matrix method to obtain the initial optimized values of the distance correction coefficient and the angle correction coefficient. Then the initial optimized values are substituted into the total error function to calculate the current total error. If the total error is greater than the preset error threshold, the error threshold is set to 5% of the standard target theoretical reflection intensity, the distance correction coefficient and the angle correction coefficient are adjusted by the initial optimized values as the starting point and by a preset step (the step value is 0.01), the total error is recalculated and the above process of calculating the partial derivative and solving the equation set is repeated until the total error is less than or equal to the error threshold. At this time, the distance correction coefficient and the angle correction coefficient obtained are the final optimized results. In the model implementation process, for each three-dimensional point in the preliminary purified data set, the actual measured distance and the incident angle are extracted, substituted into the intensity normalization model, and the corrected intensity value is calculated, that is, the corrected intensity value = original intensity value x distance correction coefficient x (square of actual measured distance reference value / square of actual measured distance) x angle correction coefficient x cos (incident angle), wherein the actual measured distance reference value is set to 100 meters. The intensity value adjusted by the model eliminates the distortion caused by the difference in actual measured distance and the change in incident angle, and generates an intensity normalized data set.

[0077] The intensity normalized data set is subjected to coordinate correction, the geometric deviation caused by the attitude change of the unmanned aerial vehicle platform is compensated based on the sensor pose parameters and the inertial measurement unit data, and a geometric correction data set is generated. Specifically, when the intensity normalized data set is subjected to coordinate correction, the three-dimensional position parameters of the laser radar sensor at the data collection time are obtained from the global navigation satellite system carried by the unmanned aerial vehicle, the position parameter accuracy is centimeter level, the attitude parameters such as pitch angle, roll angle and yaw angle of the sensor at the same time are obtained from the inertial measurement unit, the attitude parameter accuracy is 0.1 degree level, the coordinate conversion matrix is constructed based on the translation transformation of the position parameters and the rotation matrix of the attitude parameters, the original coordinates of each three-dimensional point in the intensity normalized data set are converted and operated by using the matrix, the geometric deviation caused by the attitude change in the flight process of the unmanned aerial vehicle is compensated, and a geometric correction data set is generated.

[0078] The geometric correction data set is sampled and optimized, the point density is adaptively adjusted based on the slope terrain features, the redundant data points are removed, and the optimized laser radar observation data set is generated, specifically including: when the geometric correction data set is sampled and optimized, the three-dimensional point cloud data in the geometric correction data set is partitioned and processed according to a 10m*10m grid, the slope value of the three-dimensional points in each partition is calculated through the elevation difference and the horizontal distance of adjacent points, the curvature value of each partition is calculated, the target point density of different partitions is set according to the slope value and the curvature value, wherein the target point density of the partition with a slope value greater than 30 degrees or a curvature value greater than 0.05 per meter is set to 20 points per square meter to 30 points per square meter, the target point density of the partition with a slope value of 15 degrees to 30 degrees and a curvature value of 0.02 per meter to 0.05 per meter is set to 10 points per square meter to 20 points per square meter, and the target point density of the partition with a slope value less than 15 degrees and a curvature value less than 0.02 per meter is set to 5 points per square meter to 10 points per square meter, the three-dimensional points in each partition are uniformly sampled and screened according to the target point density, the key feature points are retained and the redundant data points exceeding the target density are removed, and the optimized laser radar observation data set is generated.

[0079] In the embodiment of the application, the discrete noise is identified by statistical characteristics, the value range is adapted to different slope point cloud densities, the terrain feature points can be effectively retained while reducing noise interference. The step corrects the intensity distortion through physical characteristic modeling, the value range covers the conventional monitoring height and angle of the unmanned aerial vehicle, and the consistency of the intensity features of different regional ground objects can be improved.

[0080] As shown in FIG. 1, Figure 2 In another preferred embodiment of the application, stable terrain feature points in the observation data set are extracted and matched, a feature point matching algorithm is used to calculate the spatial transformation parameters, the optimized laser radar observation data is accurately registered, the conversion and fusion of all data to a unified spatial reference are completed, and a slope three-dimensional observation data set is generated, including:

[0081] In the observation data, the stable terrain feature points are extracted, the surface protrusions, the rock layer junction lines and the terrain mutation points are selected as the high reliability feature points based on the curvature and the normal vector change rate indexes, and an initial feature point set is generated. Specifically, when the stable terrain feature points are extracted in the observation data, for each three-dimensional point in the optimized laser radar observation data set, a neighborhood point set with the point as the center and a radius of 0.5 meters to 2 meters is selected, a local quadratic surface is fitted through the neighborhood point set, the curvature value of the point is calculated, the curvature value is calculated by using the average value of the principal curvatures of the surface, and the value range is 0 to 0.5 per meter; and the normal vector angle of the point and all points in the neighborhood point set is calculated, the normal vector change rate is obtained by the ratio of the angle difference value and the distance between points, and the normal vector change rate value range is 0 to 30 degrees per meter; the curvature threshold is set to 0.1 per meter, the normal vector change rate threshold is set to 15 degrees per meter, the points with the curvature value greater than or equal to the curvature threshold or the normal vector change rate greater than or equal to the normal vector change rate threshold are selected, and are marked as the surface protrusions, the rock layer junction lines or the terrain mutation points, and the initial feature point set is generated.

[0082] When the least square method is used to fit the surface of the sample points, a quadratic surface model to be fitted is determined, the model includes a quadratic term of x coordinate, a quadratic term of y coordinate, a cross term of x and y, a first order term of x coordinate, a first order term of y coordinate and a constant term, and the model is used to describe the terrain surface shape of the local area where the sample points are located. When the sample points are selected, a neighborhood with a radius range of 0.5 meters to 2 meters is set with the target three-dimensional point as the center, all three-dimensional points in the neighborhood are extracted as the fitting sample points, the number of sample points is controlled to be 50 to 200, and enough terrain feature information is ensured to be covered. When the error function is constructed, the square of the difference between the actual z coordinate value of each sample point and the predicted z coordinate value calculated by the surface model is used as the single sample error item, the total error function is obtained by accumulating the single sample error items of all sample points, and the total error function reflects the deviation degree of the model prediction value and the actual terrain.

[0083] The partial derivatives of the total error function with respect to each coefficient in the model are calculated, and linear equations are established based on the extreme condition that the partial derivatives are zero. The unknowns in the equations are the coefficients of the surface model, the coefficient matrix of the equations is composed of the x coordinate values, y coordinate values and constant term of the sample points, and the constant term vector is composed of the product sum of the actual z coordinate values of the sample points and the coordinate combination values. The linear equations are solved by matrix inversion operation to obtain the initial values of the coefficients of the surface model, and the initial coefficients are substituted into the surface model to calculate the predicted z coordinate values of each sample point. The residuals between the actual z coordinate values and the predicted z coordinate values are calculated, the standard deviation is calculated based on the residuals, and if the standard deviation is greater than a preset threshold (the threshold is set to 0.05 meters), the abnormal sample points with absolute residual values greater than twice the standard deviation are removed, and the remaining sample points are selected to repeat the above fitting process until the residual standard deviation is less than or equal to the threshold or the iteration number reaches a preset upper limit (the upper limit of the iteration number is set to 5 times), at which time the obtained coefficients are the final coefficients of the final surface model.

[0084] The least squares method is used to construct the total error function E, which is the sum of the error squares of the actual z coordinate values and the predicted z coordinate values of all sample points, and its calculation formula is where n is the total number of sample points, is the actual z coordinate value of the i-th sample point, is the predicted coordinate value of the i-th sample point; substituting into the above formula, the total error function can be expanded as E= where , is the plane coordinate value of the i-th sample point, , , , , , are the coefficients to be solved of the quadratic polynomial model; after expansion, it is further expressed as E= + ; where is the horizontal coordinate linear contribution term of the i-th sample point; is the vertical coordinate linear contribution term of the i-th sample point; is the horizontal coordinate nonlinear contribution term of the i-th sample point; is the horizontal and vertical coordinate coupling contribution term of the i-th sample point; is the vertical coordinate nonlinear contribution term of the i-th sample point.

[0085] , respectively represent the x coordinate value and y coordinate value (horizontal direction coordinate) of the i-th sample point in the three-dimensional space, which are used to describe the plane 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. Represents the first-order term of the x-coordinate in the quadratic surface model ( 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 sample point selection needs to determine the range in combination with the complexity of the slope terrain. A circular neighborhood with a radius value range of 0.5 meters to 2 meters is set with the target three-dimensional point as the center. In the rock area, the neighborhood radius is 0.5 meters to 1.2 meters due to the relatively regular terrain. In the soil area, the neighborhood radius is 1.2 meters to 2 meters due to the loose soil and the complex terrain caused by the presence of stones. It is ensured that the neighborhood contains a sufficient number of sample points to reflect the local terrain features. The number of sample points is controlled to be 50 to 200. The number of sample points in a single neighborhood in the rock area is 50 to 120, and the number of sample points in the soil area is 120 to 200. If the number of sample points in the neighborhood is less than 50, the neighborhood radius is expanded until the number of sample points meets the standard. If the number of sample points exceeds 200, some redundant points are removed by uniform sampling (every 2 points, keep 1), to avoid excessive data increasing the calculation load while ensuring the representativeness of the sample points.

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

[0089] The initial feature point set is constructed with a feature descriptor, and a feature description vector is generated based on the local surface geometric properties to obtain a standardized feature description set. The standardized feature description set is matched for similarity, and the mis-matched point pairs are removed through bidirectional consistency inspection to generate a reliable feature point correspondence set, which specifically includes: when constructing the initial feature point set with a feature descriptor, for each feature point, a local neighborhood with a radius of 1 meter to 3 meters is selected with the point as the center. The average curvature, normal vector direction, elevation standard deviation and intensity mean of all points in the neighborhood are calculated. The geometric properties are arranged in a predetermined order to form a basic feature vector. The basic feature vector is standardized to convert each attribute value to the interval range of 0 to 1. The dimension difference of different attributes is eliminated by linear scaling to form a standardized feature description vector with a length of 128 dimensions. The description vectors of all feature points form the standardized feature description set.

[0090] Based on the set of reliable feature point correspondences, the final spatial transformation parameters are calculated to generate a set of coordinate transformation parameters; the set of coordinate transformation parameters is used to rigidly transform all lidar observation data to be uniformly converted to the reference coordinate system to generate a set of uniformly spatially referenced data, specifically including: when calculating the final spatial transformation parameters based on the set of reliable feature point correspondences, the point pairs in the set of reliable feature point correspondences are taken as input to construct a spatial transformation model composed of a rotation matrix and a translation vector, and the initial transformation parameters are calculated by the three-point method. An iterative optimization method is used to minimize the distance error of the corresponding points after transformation, and the Euclidean distance of all corresponding points after transformation is calculated each time. Points with a distance greater than the sum of the average distance and twice the standard deviation are marked as outliers and temporarily excluded. Based on the remaining inliers, the transformation parameters are recalculated, and the iteration process is repeated until the variation of the transformation parameters in the last two iterations is less than 0.01 radian and 0.05 meters. At this time, the rotation matrix and the translation vector form the set of coordinate transformation parameters; when using the set of coordinate transformation parameters to rigidly transform all lidar observation data, the three-dimensional point coordinates in each period of data are substituted into the rotation matrix and the translation vector in the set of coordinate transformation parameters. The angle deviation of different periods of data is eliminated by rotation operation, and the spatial position reference of the data is unified by translation operation. The reference coordinate system is selected as the coordinate system of the first period of observation data, and the coordinates of all data points after transformation are converted to the reference coordinate system to generate a set of uniformly spatially referenced data.

[0091] 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, which contains two types of core parameters: rotation and translation. The rotation parameters describe the spatial rotation relationship through three Euler angles (pitch angle, roll angle, yaw angle), and the translation parameters describe the spatial position adjustment relationship through three-dimensional coordinate offset (Δx, Δy, Δz). The mathematical expression of the model is that the three-dimensional coordinates in the target coordinate system are equal to the coordinates in the original coordinate system after rotation matrix conversion and superposition of translation vector, so as to realize the spatial position alignment of lidar data in different periods.

[0092] The model training stage takes a reliable feature point correspondence set as an input sample, and 30 to 100 pairs of feature points are randomly selected from the correspondence set as training sample pairs. The sample pairs are evenly distributed in the rock area, soil area and boundary area of the slope to avoid being concentrated in a single land element. When establishing the error function, the Euclidean distance square of each pair of feature points calculated by the transformation model is taken as the single sample error term, and the error terms of all sample pairs are accumulated to obtain the total error function, which reflects the overall fitting effect of the transformation model. The partial derivatives of the total error function with respect to the rotation parameters and the translation parameters are calculated, and a nonlinear equation set is established according to the extreme condition that the partial derivative is zero. The equation set is solved by 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 correction amount is calculated by substituting the equation set, and the total error is recalculated after updating the rotation matrix and the translation vector. The iteration process is repeated until the total error difference between the previous and the current iteration is less than 0.01 meters or the iteration number reaches 20, and the preliminary optimized transformation parameters are obtained. The preliminary optimized parameters are purified to remove outliers, the residual Euclidean distance of all feature point pairs after transformation is calculated, the average value and the standard deviation of the residual are calculated, and the feature point pairs with a residual greater than the average value plus twice the standard deviation are marked as abnormal matching pairs and removed. The remaining feature point pairs are retained to perform the above training process again to generate the final optimized spatial transformation parameter set, wherein the rotation angle accuracy is controlled within 0.01 degrees, and the translation vector accuracy is controlled within 0.05 meters. In the model implementation process, the optimized rotation matrix and translation vector are applied to all periods of laser radar observation data, and the original coordinates of each three-dimensional point are rotated and translated to convert them to the reference coordinate system. After the conversion is completed, 50 to 100 feature points in different regions of the slope are randomly selected to verify the registration accuracy, the average distance error of the same named points after registration is calculated to ensure that the error is less than 0.1 meters. If the error exceeds the threshold, the training sample pairs are reselected for model training and optimization until the accuracy requirements are met.

[0093] The data sets unified in space are fused to eliminate redundant points in the overlapping area and generate a complete three-dimensional observation data set of the slope, which specifically includes: when the data sets unified in space are fused, the data sets are partitioned according to a 5m x 5m grid, the number of three-dimensional points in each grid is counted, and when the number of points in a grid exceeds a preset density threshold (the threshold is set to 50 points per grid), the point cloud is retained by using the uniform sampling method, i.e. every 0.5 meters is retained according to the spatial distribution of the points in the grid, and the redundant points exceeding the density threshold are removed. When the number of points in a grid is less than the density threshold, all points are retained to ensure the integrity of the terrain. The fused data covers the entire slope area without overlapping and redundancy, generating a complete three-dimensional observation data set of the slope.

[0094] In the embodiment of the present application, the standardized feature description set is constructed by local surface geometric properties, and the false matching point pairs are removed by combining with the bidirectional consistency check, so that the matching error caused by similar features in complex terrain can be reduced. In the observation data, the stable terrain feature points are selected based on the curvature and normal vector change rate indexes, so that the feature points with geometric recognition degree such as surface protrusions, rock layer junction lines and terrain mutation points can be accurately identified, and the unstable points susceptible to environmental interference can be removed.

[0095] In a preferred embodiment of the present application, based on the three-dimensional observation data set of the slope, two stable reference bases are preset, a directional reference frame is established, and adaptive correction factors are generated according to the topological relationship and structural features of each region in the frame, including:

[0096] The first space reference point located in the stable bedrock area at the top of the slope and the second space reference point located in the stable bedrock area at the toe of the slope are selected from the three-dimensional observation data set of the slope, and an initial reference axis is generated, specifically including:

[0097] The initial reference axis is taken as the reference direction, a regional coordinate system with the first space reference point as the origin and the reference line direction as the main axis is constructed, and a directional reference frame is generated, specifically including: the mountainous highway specific section slope carries out comprehensive pre-geological investigation, and the bedrock distribution of the top of the slope and the toe area is clear, the area without bedding sliding hidden danger at the top of the slope and the bedrock integrity meeting the requirements is determined as the candidate area of the first space reference point, and the area without loose silty clay accumulation and the bedrock without weathering risk at the toe of the slope is determined as the candidate area of the second space reference point.

[0098] The stereoscopic observation data set of the slope is acquired by combining a UAV oblique photogrammetry technique with a ground three-dimensional laser scanning technique, wherein the UAV oblique photogrammetry technique is used to acquire three-dimensional terrain data of the whole slope, and the ground three-dimensional laser scanning technique is used to acquire high-precision point cloud data of the candidate area at the top and the toe of the slope. Then, the point cloud data of the candidate area is denoised to remove invalid point cloud data caused by tree obstruction and dust interference, and to retain effective point cloud data reflecting the actual morphology of the bedrock. Three non-collinear points are selected from the effective point cloud data of the candidate area at the top of the slope after processing, and the coordinates of the three points are measured in the field by a total station, and the average value of the coordinates of the three points is calculated as the three-dimensional coordinates of a first spatial reference point, so as to ensure that the first spatial reference point is located in the stable bedrock area at the top of the slope and can represent the spatial position reference of the area. Meanwhile, three non-collinear points are selected from the effective point cloud data of the candidate area at the toe of the slope after processing, and the coordinates of the three points are also measured in the field by a total station, and the average value of the coordinates of the three points is calculated as the three-dimensional coordinates of a second spatial reference point, so as to ensure that the second spatial reference point is located in the stable bedrock area at the toe of the slope and can represent the spatial position reference of the area; the three-dimensional coordinate calculation software is used to connect the three-dimensional coordinates of the first spatial reference point and the three-dimensional coordinates of the second spatial reference point, so as to form a straight line connecting the two reference points, and the straight line is the initial reference axis.

[0099] Based on the directional reference frame, the monitoring area is divided into a plurality of independent analysis units according to the geological structure characteristics and the topographic change trend of the slope, and a unit division result is generated, specifically including: the first spatial reference point is set as the origin of a regional coordinate system, and the initial reference axis is set in the direction from the first spatial reference point to the second spatial reference point as the main axis direction of the regional coordinate system, which is consistent with the length direction of the slope to cover the whole 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 main axis direction is determined in combination with the topographic slope of the slope, which is parallel to the transverse section direction of the slope and is used to reflect the spatial position in the width direction of the slope.

[0100] Determine the vertical coordinate axis direction of the regional coordinate system, which is perpendicular to the local horizontal plane, and obtain the reference height of the local horizontal plane by leveling to ensure that the vertical coordinate axis direction accurately reflects the height direction of the slope; calibrate the three coordinate axes of the regional coordinate system to determine the correspondence between the unit length of each coordinate axis and the actual physical length, and set the unit length to a length unit that can accurately reflect the small deformation of the slope to meet the monitoring accuracy requirements; then convert all three-dimensional data in the slope stereoscopic observation data set to the regional coordinate system, convert the coordinate values of the original observation data to the coordinate values in the regional coordinate system through the coordinate conversion algorithm, and ensure that all monitoring data are processed in the same coordinate system; organize the origin position coordinate axis direction unit length and data conversion rules of the regional coordinate system to form a directional reference framework, which can accurately correspond to the actual spatial position and topographic features of the slope, and provide a unified spatial reference basis for subsequent analysis unit division.

[0101] Extract the geological structure feature data of the slope in the directional reference framework, which includes the distribution range of the upper moderately weathered sandstone interbedded with mudstone, the angle between the rock layer tendency and the slope direction, and the distribution range of the lower Quaternary residual slope deposit silty clay, the soil loose degree and the block stone distribution, and extract the topographic change trend data of the slope, including the slope gradient change of different regions of the slope, the position distribution of the slope shoulder, waist and foot, and the topographic relief characteristics; secondly, determine the division boundary according to the extracted geological structure feature data, set the boundary between the upper moderately weathered sandstone interbedded with mudstone and the lower silty clay as the longitudinal division boundary, and set the boundary between the area with bedding sliding hidden danger and the area without bedding sliding hidden danger as the transverse division boundary.

[0102] Determine the division interval in combination with the topographic change trend data, use smaller division interval in the silty clay area with steep slope and easy collapse and the sandstone interbedded with mudstone area with prominent bedding sliding hidden danger to improve the density of unit division, and use larger division interval in the area with gentle slope and stable geological conditions to reduce redundant units; according to the determined division boundary and division interval, divide the entire slope monitoring area into multiple independent analysis units in the directional reference framework, and the spatial range of each independent analysis unit can completely contain a specific geological structure feature and topographic change feature, ensuring the consistency of the geological and topographic conditions in each unit; number each independent analysis unit, record the four vertex coordinates of each unit in the directional reference framework, and describe the geological structure type and topographic feature corresponding to the unit, and organize the information to form the unit division result.

[0103] extracting terrain structure attributes of each analysis unit to generate a unit terrain feature description set; based on the relative spatial position relationship between each analysis unit and two spatial reference points, calculating the orientation weight and distance weight of each unit relative to the reference line to generate a unit weight distribution set, specifically including: for each independent analysis unit in the unit division result, extracting the terrain structure attributes of the unit through a three-dimensional terrain analysis software, the extracted terrain structure attributes including the average slope, maximum slope, minimum slope, average slope direction, slope direction variation range, ground relief of the unit, and the stratum type corresponding to the unit, if the unit belongs to the lower silty clay stratum, the compactness grade and block stone content ratio of the silty clay in the unit also need to be extracted, these attribute data are obtained by statistical analysis on the point cloud data of the corresponding unit in the slope three-dimensional observation data set. Then the terrain structure attributes of each independent analysis unit are arranged in a predetermined format, each unit corresponds to an attribute record, and the record content includes unit number, average slope, maximum slope, minimum slope, average slope direction, slope direction variation range, ground relief, stratum type compactness grade and block stone content ratio. The attribute records of all units are summarized to form a unit terrain feature description set.

[0104] The orientation weight of each independent analysis unit relative to two spatial reference points is calculated. First, the geometric center point coordinates of each unit are determined in the directional reference frame, and the included angle between the geometric center point and the initial reference axis is calculated, which is the azimuth angle of the unit. According to the risk characteristics of the slope, the azimuth angle weight mapping rule is set, wherein the units located in the azimuth angle range of the bedding sliding hidden danger section and the units located in the azimuth angle range of the silty clay easy collapse section are given higher azimuth weight values, and the units located in the azimuth angle range of the stable geological condition area are given lower azimuth weight values. The azimuth weight of each unit is calculated according to the mapping rule; the distance weight of each independent analysis unit relative to two spatial reference points is calculated, the straight line distance from the geometric center point of each unit to the first spatial reference point and to the second spatial reference point is calculated respectively, and the distance weight calculation rule is set according to the correlation between distance and reference point stability. The units closer to the two spatial reference points are given higher distance weight values due to the influence of the stability of the reference points, and the units farther away from the two spatial reference points are given lower distance weight values. The distance weight of each unit is calculated according to the calculation rule; the unit number, azimuth weight value and distance weight value of each independent analysis unit are arranged to form a unit weight distribution set.

[0105] The unit terrain feature description set and the unit weight distribution set are fused, the terrain stability comprehensive index of each analysis unit is calculated, and the unit stability evaluation set is generated, specifically including: quantitatively processing each terrain structure attribute in the unit terrain feature description set, setting a quantitative scoring standard according to the influence degree of each attribute on the slope stability, wherein the greater the average slope is, the lower the score is, the smaller the angle between the average slope direction and the rock stratum inclination is, the lower the score is, the greater the ground relief is, the lower the score is, the lower the compactness grade of silty clay is, the higher the score is, the higher the block stone content ratio is, the higher the score is, the score of the medium weathered sandstone stratum with mudstone is higher than that of the silty clay stratum, and each terrain structure attribute of each independent analysis unit is scored according to the quantitative scoring standard to obtain the attribute score value of each unit; the azimuth weight value and the distance weight value of each independent analysis unit are extracted from the unit weight distribution set, and the azimuth weight value and the distance weight value are combined and calculated according to a preset proportion to obtain the comprehensive weight value of each unit, wherein the azimuth weight value has a higher proportion in the comprehensive weight value calculation than the distance weight value, considering that the slope bedding sliding and collapse risk are mainly related to the azimuth.

[0106] Each attribute score value of each independent analysis unit is multiplied by the comprehensive weight value of the unit to obtain a weighted score value corresponding to each attribute score, and the weighted score values corresponding to all attribute scores of each unit are summed to obtain the terrain stability comprehensive index of the unit, wherein the greater the value of the comprehensive index is, the better the terrain stability of the unit is, and the smaller the value is, the worse the terrain stability of the unit is; the unit number and the corresponding terrain stability comprehensive index of each independent analysis unit are sorted according to the unit number sequence to form the unit stability evaluation set.

[0107] 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, specifically including: determining the types of ground object classification parameters involved in slope monitoring, which include slope surface displacement monitoring parameters, slope deep displacement monitoring parameters, soil pore water pressure monitoring parameters, and ground crack development degree monitoring parameters, and the ground object classification parameters are key monitoring parameters reflecting the slope stability, and the monitoring accuracy directly affects the slope instability risk judgment result; setting an adjustment coefficient value rule according to the terrain stability comprehensive index value of each independent analysis unit in the unit stability evaluation set, for the unit with a stability level of stable, since its geological and terrain conditions are stable and the monitoring parameters are less affected by environmental interference, the corresponding ground object classification parameter adjustment coefficient is set to a smaller value to reduce unnecessary correction operation;

[0108] For the unit with the stable level of relatively stable, since it has a slight risk of stability fluctuation, the corresponding ground object classification parameter adjustment coefficient is set to a medium value to achieve moderate correction; for the unit with the stable level of unstable, since it is prone to bedding slip or collapse and the monitoring parameters are easily affected by complex environment, the corresponding ground object classification parameter adjustment coefficient is set to a larger value to ensure the accuracy of the monitoring parameters; for each independent analysis unit, the adjustment coefficient of each type of ground object classification parameter corresponding to the unit is calculated according to the stability level corresponding to the terrain stability comprehensive index of the unit and the preset adjustment coefficient value rule, so that each type of ground object classification parameter has an adjustment coefficient matched with the stability of the unit; the unit number, the name of each type of ground object classification parameter and the corresponding adjustment coefficient of each independent analysis unit are sorted and summarized to form the adaptive correction factor.

[0109] In the embodiment of the present application, through the double-reference directional and partition weighting mechanism, the unit size and weight range adapt to the terrain differences of the composite slope, which can improve the correction accuracy of the complex terrain area, and through the stability-related parameter adjustment, adaptive parameter optimization basis is provided for subsequent ground object classification processing, and the classification deviation caused by the variable terrain is reduced.

[0110] In a preferred embodiment of the present application, the ground object classification processing parameters are optimized according to the adaptive correction factor to realize accurate separation of vegetation cover, artificial structures and natural ground surface, and to generate an optimized terrain feature data set, which includes:

[0111] According to the adaptive correction factor, the parameter adjustment set of each analysis unit is generated; based on the parameter adjustment set of each unit, the classification threshold parameter in the ground object classification algorithm is adaptively adjusted according to the terrain stability characteristics of different analysis units to generate an optimized classification parameter set, which specifically includes: when the ground object classification processing parameters are optimized according to the adaptive correction factor, the adaptive correction factor is matched with each analysis unit one by one, each analysis unit corresponds to a unique parameter adjustment coefficient, the coefficient value range is 0.6 to 1.4, wherein the coefficient corresponding to the unit with high terrain stability is greater than 1.0, and the coefficient corresponding to the unit with low terrain stability is less than 1.0, and a unit parameter adjustment set containing all unit adjustment coefficients is generated; when the classification threshold parameter is adjusted based on the unit parameter adjustment set, the initial reflectance intensity threshold range for vegetation cover classification is set to 20 to 60, the initial reflectance intensity threshold range for artificial structure classification is set to 60 to 100, and the initial reflectance intensity threshold range for natural ground surface classification is set to 10 to 30, the threshold interval between vegetation and natural ground surface is expanded by 10% to 20% for high stability units through the parameter adjustment coefficient, and the threshold interval between artificial structures and natural ground surface is reduced by 5% to 15% for low stability units, and an optimized classification parameter set exclusive to each analysis unit is generated.

[0112] According to the optimized classification parameter set, the preliminary feature classification of the slope stereoscopic observation data set is performed, and an initial classification result is generated;The initial classification result is processed to generate an optimized classification result, specifically including: when the preliminary feature classification is performed according to the optimized classification parameter set, the reflection intensity value and the local slope value of each three-dimensional point in the slope stereoscopic observation data set are extracted point by point, the points with the reflection intensity value falling within the vegetation threshold range and the slope less than 25 degrees are marked as vegetation covered points, the points with the reflection intensity value falling within the artificial structure threshold range and the slope less than 15 degrees are marked as artificial structure points, and the remaining points are marked as natural ground points, and the initial classification result is generated;When the initial classification result is processed, the number of three-dimensional points in each classification region is counted, the isolated point group with the number less than 50 is determined as a classification error and is reclassified to the adjacent dominant classification region, and the edge points of the vegetation covered region and the artificial structure region are subjected to buffer analysis, the buffer distance is set to 0.5 meters, the fuzzy points in the buffer area are reclassified according to the majority principle, and the optimized classification result is generated.

[0113] Based on the optimized classification result, the spatial distribution information of the vegetation covered region and the artificial structure region is extracted, and a feature mask data set is generated;The feature mask data set is used to process the slope stereoscopic observation data set, and the data points corresponding to the vegetation and the artificial structure are removed, and a bare ground three-dimensional point set is generated, specifically including: when the feature mask data set is generated based on the optimized classification result, the three-dimensional points of the vegetation covered region and the artificial structure region are subjected to spatial rasterization processing, the grid size is set to 0.5m*0.5m, the grid containing the vegetation or artificial structure point is marked as a mask region, and the feature mask data set composed of grid coordinates and mask identification is generated;When the feature mask data set is used to process the slope stereoscopic observation data set, the mask identification of the grid where the three-dimensional point is located is checked point by point, all three-dimensional points located in the mask region are removed, and the three-dimensional points not covered by the mask are retained, and the bare ground three-dimensional point set is generated.

[0114] The bare ground three-dimensional point set is processed to generate an optimized terrain feature data set, specifically including: when the bare ground three-dimensional point set is processed, the discrete noise points are removed by adopting the neighborhood analysis with a radius of 0.3 meters, the noise point determination standard is the isolated point with less than 5 points in the neighborhood, then the point density is calculated by 10m*10m grid partition, the region with the density less than 3 points / m2 is interpolated, the region with the density higher than 20 points / m2 is uniformly down-sampled, and finally the optimized terrain feature data set with the point density controlled in the range of 5 to 15 points / m2 is generated.

[0115] In the embodiment of the application, the adaptive parameter adjustment is used to realize the accurate separation of features under complex terrain, the value range is adapted to the distribution characteristics of the slope vegetation and artificial structure, the feature classification accuracy is improved, and the interference of the vegetation and the artificial structure on the terrain deformation analysis is effectively eliminated.

[0116] In a preferred embodiment of the present application, the terrain feature dataset is spatially interpolated to construct a digital terrain surface model, comprising:

[0117] The optimized terrain feature dataset is analyzed for terrain features, terrain feature points are identified based on local curvature variation and elevation variation coefficient, 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 sampled with optimized density, the sampling density is adaptively adjusted based on terrain complexity, sparse sampling in flat terrain areas and dense sampling in complex areas to generate a structured terrain point set, specifically comprising: when analyzing the optimized terrain feature dataset, for each three-dimensional point in the dataset, a neighborhood region with a radius range of 0.5 meters to 2 meters centered on the point is selected, and 30 to 100 three-dimensional points in the neighborhood are extracted as analysis samples. The local curvature value of the point is calculated by surface fitting, the curvature value calculation range is set to 0 to 0.5 per meter, and the ratio of the standard deviation of the elevations of all points in the neighborhood to the average elevation is calculated as the elevation variation coefficient, the variation coefficient value range is 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 mutation points are identified in combination with the slope terrain strike, wherein the ridge line area needs to satisfy that the elevations of more than 3 consecutive potential feature points are higher than those of the two sides, the valley line area needs to satisfy that the elevations of more than 3 consecutive potential feature points are lower than those of the two sides, and the slope mutation point needs to satisfy that the slope difference value of adjacent points is greater than 15 degrees, and the candidate terrain feature point set is generated after screening.

[0118] When the candidate terrain feature point set is sampled with optimized density, the slope area is divided into several sub-areas according to a 10m x 10m grid, the average slope value and curvature value of each sub-area are calculated; when the average slope of a sub-area is less than 15 degrees and the average curvature is less than 0.05 per meter, it is determined as a flat terrain area, and the sampling density is set to 5 points per square meter to 10 points per square meter; when the average slope is 15 degrees to 30 degrees and the average curvature is 0.05 per meter to 0.1 per meter, it is determined as a medium complex terrain area, and the sampling density is set to 10 points per square meter to 20 points per 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 as a complex terrain area, and the sampling density is set to 20 points per square meter to 30 points per square meter; the feature points are selected in each sub-area using uniform sampling method to ensure that the sampling points uniformly cover the area and the key terrain nodes are not missed, and a structured terrain point set is generated.

[0119] The structured terrain point set is used to construct a spatial interpolation function to obtain final interpolation parameters and generate an optimized interpolation function; the optimized interpolation function is used to calculate the numerical value of the regularly distributed elevation points to generate a preliminary terrain surface model, which specifically includes: when the structured terrain point set is used to construct a spatial interpolation function, an inverse distance weighted interpolation model is selected as a basic framework, and the spatial position and elevation value of the sampling point are used as input parameters; during the construction process, the interpolation search radius is set to 5-10 meters to ensure that 8-20 structured terrain points are included around each interpolation calculation point; when determining the weight coefficient, different weight values (1-5) are selected for trial calculation, 20% of the structured terrain points of each trial result are selected as verification points, and the average absolute error of the predicted elevation and the actual elevation of the verification points is calculated; when the weight coefficient is 2-3, the average absolute error is controlled within 0.1 meters, and the weight coefficient is determined as the final parameter to generate the optimized interpolation function; when implemented, the slope region is divided into regular grids with a side length of 1-2 meters, the optimized interpolation function is applied to each grid node, the elevation value of the node is obtained by weighted calculation of the elevations of the structured terrain points within the search radius, and a preliminary terrain surface model is generated.

[0120] When the optimized interpolation function is constructed, an inverse distance weighted interpolation model is used as a basic framework, which realizes elevation prediction by calculating the distance weight of the to-be-interpolated point and the surrounding known terrain points. The closer the known point, the greater the influence weight on the interpolation result. The slope monitoring region is divided into regular calculation units according to a 5m*5m grid, and the center point of each unit is used as a to-be-interpolated point to form a regularly distributed elevation calculation grid; during the model training stage, the three-dimensional points in the structured terrain point set are used as known sample points, and a dynamic search radius is set for each to-be-interpolated point. The search radius is adaptively adjusted according to the terrain complexity of the region, and the search radius is set to 15-20 meters in a flat terrain region and 5-10 meters in a complex terrain region to ensure that 20-50 known sample points can be searched around each to-be-interpolated point.

[0121] The sample points searched are sorted by the Euclidean distance from the point to be interpolated, and the weight value of each sample point is calculated. The closer the distance, the greater the weight value, and the weight value decreases non-linearly with the increase of the distance. The weight decay coefficient is set to 1.5 to 3.0 according to the terrain roughness, and a larger decay coefficient is selected for areas with high roughness; when establishing the error evaluation mechanism, 10% of the points in the structured terrain point set are randomly selected as verification samples, which do not participate in the interpolation calculation. The elevation values of these verification samples are predicted by the interpolation function, the absolute error between the predicted elevation and the actual elevation is calculated, and the average error and the maximum error of all verification samples are counted; if the average error is greater than 0.1 meters or the maximum error is greater than 0.3 meters, adjust the search radius and the weight decay coefficient to recalculate, wherein the search radius is increased or decreased by 2 meters, and the decay coefficient is adjusted by 0.5 steps, and the 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. At this time, the search radius range and the weight decay coefficient obtained are the final interpolation parameters; in the model implementation process, the final interpolation parameters are applied to all points to be interpolated. For each calculation unit center point, the surrounding sample points are selected according to the optimized search radius, the weight value of each sample point is calculated based on the weight decay coefficient, and the elevation value of the point to be interpolated is calculated by weighted average to generate a preliminary terrain surface model. In the implementation process, each interpolation result is checked for local terrain continuity. If the elevation difference between adjacent points to be interpolated exceeds the reasonable difference value corresponding to the average slope of the region, the weight parameters are adjusted locally to ensure smooth transition of the terrain surface.

[0122] The preliminary terrain surface model is smoothed to eliminate local interpolation anomalies and generate an optimized terrain surface model; the optimized terrain surface model is subjected to integrity testing, and data gaps are filled based on the terrain continuity principle to generate a final digital terrain surface model, which includes: when smoothing the preliminary terrain surface model, a moving window of 3x3 to 5x5 is used to calculate the sliding average of the model elevation data, and the elevation value of the center pixel of the window is replaced by the average value of all pixel elevations in the window. The pixels with a difference between the window elevation value and the average value greater than 0.2 meters are marked as abnormal values, and the average value of the non-abnormal values in the window is used to replace them to eliminate abnormal protrusions or depressions caused by local interpolation fluctuations. When the integrity of the optimized terrain surface model is tested, all grid nodes are traversed, and regions with continuous 3 or more no elevation values are determined as data gap regions; the elevation values of the effective grid nodes within a 5-meter range around the gap region are used for interpolation to fill the gap. In the filling process, the elevation difference between adjacent nodes is ensured to be less than 0.3 meters to maintain terrain continuity, and a complete digital terrain surface model is finally generated.

[0123] In the embodiment of the present application, the sampling density and interpolation parameters are adaptively adjusted according to the terrain complexity, and the value range is adapted to the terrain diversity of the rock-soil composite slope, so that the elevation accuracy of the digital terrain surface model can be controlled within 0.1 meters, the key terrain features such as ridges and valleys are retained, the interpolation noise is eliminated, high-precision terrain reference data for subsequent slope displacement vector field calculation is provided, and the reliability of deformation monitoring under complex terrain conditions is improved.

[0124] In a preferred embodiment of the present application, based on the digital terrain surface model, the spatial similarity analysis is performed and the gradient descent algorithm is used for iterative optimization to calculate the slope surface displacement vector field, identify the potential sliding surface and deformation abnormal area, and generate the displacement field calculation result, including:

[0125] 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, specifically including: when the digital terrain surface model is divided into a plurality of analysis units, adaptive grid division is adopted according to the complexity of the slope terrain, 5m x 5m small analysis units are set at the junction between rock and soil areas and the slope mutation area, and 10m x 10m large analysis units are set in the relatively flat homogeneous area to ensure that the unit boundary is consistent with the trend of the terrain feature line; for each analysis unit, the corresponding regions of the baseline period digital terrain surface model and the monitoring period digital terrain surface model are selected, the elevation difference values of all corresponding points in the unit are calculated, the square sum of the elevation difference values is taken as the unit elevation residual objective function value, and the unit difference quantization data set is generated, wherein the elevation difference calculation accuracy is controlled within 0.01 meters.

[0126] Based on the unit difference quantization data set, the gradient descent algorithm is used for iterative optimization to solve, the gradient of the objective function is calculated, and the spatial position parameters of the analysis unit are gradually adjusted along the negative gradient direction to generate a unit displacement parameter set, specifically including: when the gradient descent algorithm is used for iterative optimization based on the unit difference quantization data set, the initial displacement parameter is set to zero vector, the learning rate value range is 0.01 to 0.1, and the maximum iteration number is set to 50 to 100 times; in each iteration process, the gradient of the current unit elevation residual objective function to the displacement parameter is calculated, the gradient direction is the direction in which the function value increases, the displacement parameter is adjusted in the opposite direction (negative gradient direction), and the adjustment amplitude is 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 function value difference between the two iterations is less than 0.05 meters, the iteration is stopped, otherwise the adjustment is continued until the maximum iteration number is reached, the three-dimensional displacement parameters of each analysis unit are generated, and the unit displacement parameter set is composed.

[0127] Based on the unit displacement parameter set, the displacement vector field of all analysis units is integrated, and the global displacement vector field is constructed. Through statistical consistency test processing, the optimized displacement vector field is generated, and the spatial distribution characteristics of displacement are analyzed to generate the potential sliding boundary identification result, including: when constructing the global displacement vector field based on the unit displacement parameter set, the displacement parameters of each analysis unit are assigned to the center points of the units, and the displacement vectors at the boundaries of the units are calculated by the neighborhood interpolation method to make the displacement vectors of adjacent units change continuously; when using statistical consistency test processing, the average difference between the displacement vector of each unit and the displacement vectors of 3 to 5 adjacent units is calculated. If the difference value is greater than 2 times the standard deviation of all unit displacement vectors, it is marked as an abnormal unit, and the weighted average value of the surrounding units is used to correct the displacement vector of the abnormal unit to generate the optimized displacement vector field. When analyzing the spatial distribution characteristics of displacement, the size and direction of the displacement vector of each unit are counted, and the area with a displacement vector size greater than 0.1 meters or a direction mutation angle greater than 30 degrees is marked as a potential sliding boundary to generate the potential sliding boundary identification result.

[0128] When the neighborhood interpolation method is implemented, the grid division is performed on the region to be interpolated, and the grid side length is set to 0.5 meters to 2 meters, wherein the grid side length of the rock slope region is 0.5 meters to 1 meter, and the grid side length of the soil slope region is 1 meter to 2 meters. Each grid node is used as a point to be interpolated. For each point to be interpolated, the neighborhood range is determined with the point as the center, and the neighborhood radius is adaptively adjusted according to the terrain complexity. The neighborhood radius of the flat terrain region is set to 10 meters to 15 meters, the neighborhood radius of the complex terrain region is set to 5 meters to 10 meters, and the neighborhood radius of the medium slope region is set to 7 meters to 12 meters. All sample points in the neighborhood range are extracted from the structured terrain point set, the horizontal distance between each sample point and the point to be interpolated is calculated, the effective sample points are selected by distance threshold, the sample points with a distance greater than the neighborhood radius are removed, and the number of sample points participating in interpolation is controlled to 10 to 30. The sample points after screening are subjected to outlier test, the average value and standard deviation of the sample point elevation are calculated, the sample points with an elevation value deviating from the average value by more than twice the standard deviation are marked as outliers and removed, and the remaining sample points are used as effective interpolation samples.

[0129] The weight coefficients of the effective sample points are calculated by using a distance inverse ratio weighting method, the weight of each sample point is inversely proportional to the horizontal distance of the point to the point to be interpolated, the closer the distance, the greater the weight; the weights of all sample points are normalized to make the total weight sum to 1, so as to avoid weight imbalance caused by distance difference; the elevation value of each effective sample point is multiplied by the corresponding weight and then accumulated to obtain the initial interpolation elevation of the point to be interpolated; the accuracy of the initial interpolation result is verified, 5% to 10% of the points in the structured terrain point set are randomly selected as verification points, the error between the actual elevation and the interpolation elevation of the verification points is calculated, and the average error and the maximum error are counted; if the maximum error exceeds 0.1 meters or the average error exceeds 0.05 meters, the neighborhood radius is adjusted to re-execute the above interpolation process until the error meets the requirements. The interpolation elevations of all grid nodes are integrated to generate continuous terrain surface elevation data, and the neighborhood interpolation calculation is completed.

[0130] The potential sliding boundary recognition result is analyzed to determine the spatial form and deformation evolution characteristics of the sliding surface, and a displacement field calculation result is generated, specifically including: when the potential sliding boundary recognition result is analyzed, the profile line is extracted along the potential sliding boundary trend combined with the slope geological survey data, the dip and inclination of the sliding surface are determined according to the direction and size change of the displacement vector on the profile line, wherein the inclination of the rock slope sliding surface is referenced to the rock layer inclination value, and the inclination of the soil slope sliding surface is combined with the slope angle value; when the deformation evolution characteristics are analyzed, the displacement vector field change of different monitoring periods is compared, the displacement rate is calculated, the area with a rate greater than 0.01 meters / day is marked as a deformation active area, and whether the deformation shows an accelerating trend is judged combined with the terrain stability index, and finally a displacement field calculation result containing the spatial coordinates of the sliding surface, the deformation direction, the rate and the active area distribution is generated.

[0131] In the embodiment of the present application, by adaptive unit division and gradient descent optimization, the value range adapts to the difference of composite slope terrain, the displacement calculation accuracy can be controlled within 0.05 meters, the accuracy of potential sliding surface recognition is improved, a quantitative basis is provided for slope instability early warning, and the monitoring reliability of high-risk slopes is effectively ensured.

[0132] In a preferred embodiment of the present application, the displacement field calculation result is input into a risk assessment model, and the slope stability is quantitatively evaluated through a multi-source data fusion analysis platform, including:

[0133] According to the displacement field solving 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, geological survey data and geotechnical mechanics parameters of the slope region are integrated to generate a geomechanics parameter set, specifically including: according to the displacement field solving result, the displacement vector field data and potential sliding surface information are extracted to generate the initial risk assessment input set, for a 500-meter-long slope of K12+300~K12+800 section of a mountainous highway, a displacement monitoring sampling point grid is determined at a density of 2m x 2m, a total of 12500 sampling points are arranged, based on the displacement field solving result, three-dimensional displacement coordinate data of each sampling point at a time interval of 30 days and 1 hour are obtained, displacement components of each sampling point in horizontal direction X-axis Y-axis and vertical direction Z-axis are calculated, displacement rate is calculated by means of displacement difference divided by time interval, displacement rate value range is 0.1~10mm / day, displacement direction angle is calculated by means of trigonometric function, displacement direction angle value range is 0~360°, displacement vector field data are integrated; combined with the stratum characteristics of the upper moderately weathered sandstone with mudstone and the lower quaternary residual slope deposit silty clay of the slope, a displacement mutation threshold is set, wherein the threshold value of the sandstone with mudstone region is 5mm / day, the threshold value of the silty clay region is 8mm / day, the displacement gradient threshold is set to 0.02mm / m, the region where the sampling point is located is marked as a potential sliding suspicious region if the displacement of the sampling point exceeds the corresponding regional mutation threshold or the displacement gradient exceeds 0.02mm / m.

[0134] Three-dimensional curved surface fitting is performed on the displacement data of the sampling points in each suspicious region to obtain spatial geometric parameters such as strike, dip angle and trend of the potential sliding surface, wherein the strike value range is 0~360°, the dip angle value range is 30~60°, and the trend value range is 0~360°; data are integrated in the format of sampling point number~displacement component~displacement rate~displacement direction angle~potential sliding surface number~sliding surface geometric parameter to generate the initial risk assessment input set.

[0135] Based on the initial risk assessment input set, geological survey data and geotechnical mechanics parameters of the slope region are integrated to generate a geomechanics parameter set, geological survey data of the slope are collected, including that the upper moderately weathered sandstone with mudstone stratum distribution range is K12+300~K12+800 section of the slope top to 60 meters high, lithology classification data is sandstone content 60~70% and mudstone content 30~40%, fault and joint development interval value range is 0.5~2m, stratum thickness is sandstone with mudstone thickness 50~60m and silty clay thickness 5~10m, geotechnical mechanics parameters are obtained through indoor direct shear test and triaxial compression test, wherein the cohesion of sandstone with mudstone is 25~35kPa, the internal friction angle is 28~35°, the elastic modulus is 15~25GPa, the Poisson's ratio is 0.2~0.25, and the specific gravity is 22~24kN / m 3, the cohesion of silty clay is valued at 15-20 kPa, the internal friction angle is valued at 18-22°, the elastic modulus is valued at 0.5-1.5 GPa, the Poisson's ratio is valued at 0.3-0.35, and the unit weight is valued at 19-21 kN / m³; according to the initial input of the strike and dip of the potential sliding surface, the range of strata through which each sliding surface passes is determined, such as a bedding sliding surface mainly passing through a sandstone interbedded with mudstone stratum and a collapse sliding surface mainly passing through a silty clay stratum; the extracted geological survey data of the corresponding stratum and the geotechnical mechanics parameters are associated according to the potential sliding surface number, and a total of 10 potential sliding surfaces are divided; the abnormal values in the associated data are detected by using the Grubbs test method with a level of 0.05, and the abnormal values detected are corrected by using the inverse distance weighted interpolation method, and the correction radius is set to 5 meters; the data is arranged according to the potential sliding surface coding classification mode, and a set of geomechanics parameters is generated.

[0136] The geomechanics parameter set and the real-time monitored hydro-meteorological data are fused to generate a multi-source fusion database; based on the multi-source fusion database, the safety factor of each potential sliding surface is calculated to generate a preliminary stability evaluation result, which specifically includes: fusing the geomechanics parameter set and the real-time monitored hydro-meteorological data to generate a multi-source fusion database, determining the hydro-meteorological monitoring indicators for the slope, including the groundwater level elevation, the pore water pressure, the rainfall, the atmospheric temperature and the relative humidity, wherein the groundwater level elevation is valued in the range of 100-185 meters with the regional datum sea level as the 0 point, the pore water pressure is valued in the range of 0-50 kPa, the rainfall is valued in the range of 0-200 mm / day, the atmospheric temperature is valued in the range of -5-35°C, and the relative humidity is valued in the range of 30-100%; three sets of groundwater level meters and pore water pressure meters are arranged at the top, middle and bottom of the slope, respectively, two sets of rain gauges and thermohygrographs are arranged at the top of the slope, and real-time data is obtained at a frequency of collecting the groundwater level and pore water pressure every 1 hour, the rainfall every 5 minutes, the temperature and humidity every 5 minutes, and the temperature and humidity every 30 minutes.

[0137] The real-time data is preprocessed, missing values are filled by using linear interpolation method, the length of single missing is not more than 2 hours, and the sliding average method with a window size of 5 data points is used to filter noise; according to the spatial coordinates of each potential sliding surface in the geomechanics parameter set, the monitoring area with a radius of 50 meters corresponding to each sliding surface is determined, and the geomechanics parameters and the hydro-meteorological data of the corresponding area are spatio-temporally aligned according to the second-level time stamp; a data association index is established with the potential sliding surface number-monitoring time-stratum number as the index keyword, the data units are unified into international standard units such as pressure unit kPa and length unit m, and the storage format is unified into JSON format; the standardized data is stored in a MySQL relational database to generate a multi-source fusion database.

[0138] The safety factor of each potential sliding surface is calculated based on the multi-source fusion database to generate a preliminary stability evaluation result. The Bishop method is selected as the safety factor calculation method. According to the strike, dip and trend of each potential sliding surface in the multi-source fusion database, the cohesion, internal friction angle, and pore water pressure data of the corresponding stratum in the database are substituted into the calculation model. The potential sliding surface is divided into several calculation blocks along the strike direction at a length of 5-8 meters. A single sliding surface is divided into 8-15 blocks. The weight of each block is calculated by multiplying the block volume by the specific gravity. The seepage force generated by the pore water pressure is calculated by multiplying the pore water pressure by the block area. The normal stress and shear stress on the sliding surface are calculated according to the geometric size and stress state of the block. The normal stress is the normal component of the weight minus the pore water pressure. The shear stress is the shear component of the weight plus the shear component of the seepage force. According to the Mohr-Coulomb strength criterion, the sliding resistance is calculated by multiplying the cohesion by the normal stress plus the product of the normal stress and the tangent of the internal friction angle. The sliding force is calculated by the shear stress. The safety factor of the entire sliding surface is solved by the iterative method. The iteration number is not more than 20 times, and the convergence error is less than 0.001. The safety factor evaluation threshold is set. When the safety factor is greater than or equal to 1.25, the state is stable. When the safety factor is between 1.05 and 1.25, the state is basically stable. When the safety factor is between 1.0 and 1.05, the state is unstable. When the safety factor is less than 1.0, the state is unstable. According to the comparison between the safety factor calculation value and the threshold, the stability state is determined when the safety factor calculation value ranges from 0.8 to 1.5. The data is arranged in the format of the potential sliding surface number to generate the preliminary stability evaluation result.

[0139] Based on the preliminary stability evaluation result, the stress and strain analysis result is generated. The stress and strain analysis result and the preliminary stability evaluation result are fused to calculate the comprehensive stability index by a weighted fusion algorithm to generate a slope stability grade evaluation set. Specifically, based on the preliminary stability evaluation result, the stress and strain analysis result is generated. The stress and strain analysis is carried out by combining field monitoring and theoretical calculation. According to the stratum distribution and lithology data and potential sliding surface information of the slope in the multi-source fusion database, stress sensors are arranged at each potential sliding surface area and stratum boundary of the slope. The sensor arrangement density is 1 per 100 square meters, and a total of 50 sensors are arranged. The tetrahedral element mesh division logic is used to determine the sensor monitoring area. The monitoring area size of the potential sliding surface and the stratum boundary area is set to 3m×3m×3m, and the monitoring area size of other areas is set to 5m×5m×5m. The elastic modulus, Poisson's ratio, and specific gravity of each stratum in the database are used as the theoretical calculation basis parameters. For the sliding surface area preliminarily evaluated as unstable and unstable, the friction characteristic parameter of the monitoring area is determined according to the tangent of the internal friction angle of the corresponding stratum. The friction characteristic parameter value ranges from 0.33 to 0.70. At the same time, the normal stress reference value is set to 5×10 5 ~1×10 6 N / m³.

[0140] In the monitoring process, the boundary stress state of the position of each sensor is recorded, the displacement of the bottom position is limited in the X-axis Y-axis Z-axis direction, the displacement of the side position is limited in the X-axis Y-axis direction, the top position has no displacement limitation, and the uniform load size of each monitoring area is determined according to the real-time pore water pressure data in the database; the stress data of the slope area is collected in real time through the sensor, the strain data is derived by combining theoretical calculation, the distribution of the maximum principal stress, the minimum principal stress and the shear stress, and the distribution of the axial strain, the shear strain and the volume strain are obtained, wherein the maximum principal stress value range is 0-200kPa, the minimum principal stress value range is -50-100kPa, the shear stress value range is 0-80kPa, the axial strain value range is -5x10 -4 -4 , the shear strain value range is 0-3x10 -4 , and the volume strain value range is -3x10 -4 -3x10 -4 ; the maximum shear stress value, the stress concentration coefficient and the maximum shear strain value of each potential sliding surface area are extracted, wherein the stress concentration coefficient is the ratio of the maximum shear stress to the average shear stress, and the value range is 1.2-2.5, the stress and strain distribution cloud diagram is drawn, the characteristic parameters and the distribution cloud diagram are integrated, and the stress and strain analysis result is generated.

[0141] The stress and strain analysis result and the preliminary stability evaluation result are fused, a comprehensive stability index is calculated by a weighted fusion algorithm to generate a slope stability grade evaluation set, four evaluation indexes are determined, which are the safety coefficient A in the preliminary evaluation result, the maximum shear stress value B in the stress and strain analysis result, the stress concentration coefficient C, and the maximum shear strain value D; the weights are determined by using the analytic hierarchy process, a judgment matrix is constructed, the characteristic vector is calculated by the sum product method, and the weights are determined after consistency test, the consistency test requires that the consistency ratio CR is less than 0.1, and the weights are A=0.5, B=0.2, C=0.2 and D=0.1; the indexes are standardized, the safety coefficient is positively standardized, the standardized value is calculated in the form of the difference between the actual value and 1.0 divided by the difference between 1.5 and 1.0, the value range is 0-1, the maximum shear stress value is negatively standardized, the standardized value is calculated in the form of the difference between 80 and the actual value divided by the difference between 80 and 0, the value range is 0-1, the stress concentration coefficient is negatively standardized, the standardized value is calculated in the form of the difference between 2.5 and the actual value divided by the difference between 2.5 and 1.2, the value range is 0-1, and the maximum shear strain value is negatively standardized, the standardized value is calculated in the form of the difference between 3x10 -4 and the actual value divided by the difference between 3x10 -4 ​Subtract the difference value of 0, the value range is 0-1; according to the formula of comprehensive stability index equal to A*0.5+B*0.2+C*0.2+D*0.1, the index value of each potential sliding surface is calculated, and the comprehensive stability index value range is 0-1; set the grade range, the comprehensive stability index is greater than or equal to 0.8, which is the first stable grade; 0.6 is less than or equal to the comprehensive stability index and less than 0.8, which is the second basic stable grade; 0.4 is less than or equal to the comprehensive stability index and less than 0.6, which is the third sub-stable grade; the comprehensive stability index is less than 0.4, which is the fourth unstable grade; according to the comparison of the calculated value and the grade range, the stability grade of each sliding surface is determined, the data is arranged according to the format of potential sliding surface number, and the slope stability grade evaluation set is generated.

[0142] Based on the slope stability grade evaluation set, the slope risk area is divided and the early warning grade is determined, and the slope risk zoning and early warning scheme is generated; according to the slope risk zoning and early warning scheme, the slope stability quantitative evaluation report is generated, which specifically includes: based on the slope stability grade evaluation set, the slope risk area is divided and the early warning grade is determined to generate the slope risk zoning and early warning scheme, according to the actual length of the K12+300-K12+800 section of the slope, it is divided into 10 evaluation units according to 50 meters as a unit, each unit corresponds to K12+300-K12+350 to K12+750-K12+800, and each unit contains 1-2 potential sliding surfaces; determine the dominant stability grade of the evaluation unit, when there is a fourth unstable grade sliding surface in the unit, the dominant grade is the fourth grade, when there is no fourth grade but there is a third grade, the dominant grade is the third grade, and so on; according to the dominant grade, the risk area is divided, the first stable corresponds to the low risk area, the second basic stable corresponds to the lower risk area, the third sub-stable corresponds to the medium risk area, and the fourth unstable corresponds to the high risk area.

[0143] Combined with the displacement rate data of the displacement vector field in the multi-source fusion database, the warning level standard was set, the low-risk area and the displacement rate <0.5 mm / day were blue warning, the lower-risk area or 0.5 mm / day≤displacement rate<1.0 mm / day was yellow warning, the medium-risk area or 1.0 mm / day≤displacement rate<2.0 mm / day was orange warning, and the high-risk area or displacement rate≥2.0 mm / day was red warning; the warning level was determined according to the dominant level and the displacement rate of each evaluation unit, the displacement rate value range was 0.1-5.0 mm / day, the risk area spatial distribution map was drawn by using the GIS software, and the range of each area was marked, such as the K12+400-K12+450 section was a medium-risk area, the dominant level, the warning level and the determination basis were 1.03, 1.2 mm / day, respectively, and the slope risk zoning and warning scheme were integrated. The beneficial effects of the process were that the evaluation unit was divided according to the actual road section to ensure that the zoning was in line with the actual engineering, the double determination of the dominant level and the displacement rate improved the warning accuracy, the GIS visualization was convenient for intuitive understanding of the risk distribution, the problems of traditional risk zoning, fuzzy warning and untimely warning were solved, and clear targeting was provided for disaster prevention and control of the high-risk slope top with houses and slope foot adjacent to the expressway.

[0144] According to the slope risk zoning and warning scheme, a slope stability quantitative evaluation report was generated, the report structure framework was determined, and six parts including project overview, data source and processing, detailed process of each evaluation step, evaluation result, risk prevention and control suggestion, and conclusion were included; in the project overview part, the geographical location of the slope was K12+300-K12+800 section of a mountainous highway, the scale was 500 meters in length, the maximum slope height was 85 meters, the average slope was 55°, the stratum characteristics were upper moderately weathered sandstone with mudstone and lower quaternary residual slope silty clay, the surrounding environment was the slope top with residential houses and the slope foot adjacent to the main line of the expressway, and the risk characteristics were bedding sliding and rain season collapse hidden danger; in the data source and processing part, it was explained that the displacement field calculation data came from the slope monitoring system with a 2m×2m sampling grid, the geological exploration data came from 8 drillings, the hydrological and meteorological data came from the monitoring equipment arranged on site including 3 sets of underground water level meters, and the pretreatment method of each data was linear interpolation and moving average.

[0145] In the detailed process part of each evaluation step, the generation process of the stress-strain analysis result grade evaluation set risk zoning scheme of the initial input set geomechanics parameter set multi-source fusion database preliminary evaluation result is described in turn, the method used in each step is clarified, such as Kriging interpolation, Bishop method, analytic hierarchy process, the parameter value is clarified, such as cohesion 25-35 kPa, safety factor threshold 1.05, and the calculation details are clarified, such as strip and block division 5-8 meters; in the evaluation result part, the safety factor of each potential sliding surface is presented, such as sliding surface 1 is 1.12 and sliding surface 5 is 0.98, the comprehensive stability index is presented, such as sliding surface 1 is 0.75 and sliding surface 5 is 0.38, the stability grade, and the risk area distribution is presented, such as K12+400-K12+500 section is a high-risk area, the early warning grade is presented, such as the section is red early warning, and the corresponding GIS distribution map is presented; in the risk prevention and control suggestion part, the anchor rod support anti-slide pile setting and drainage system optimization are proposed for the red early warning area, such as K12+400-K12+500 section, wherein the anchor rod length is 8-10 meters, the interval is 2m*2m, the pile length is 15-20 meters, the pile diameter is 1.2 meters, and the interval is 5 meters, and the drainage system optimization includes adding a water intercepting ditch and a blind ditch; the monitoring frequency is suggested to be increased to 15 minutes / second for the orange early warning area.

[0146] The suggestions of regular monitoring 1 hour / second and regular inspection 1 time per week are respectively proposed for the yellow and blue early warning areas; in the conclusion part, the overall stability state of the slope is summarized, such as 30% area is high risk and 20% is medium risk, the key risk area is clarified, such as K12+400-K12+500 section, and the priority prevention and control area is clarified, the applicable conditions of the evaluation result are clarified, such as no significant change in stratum in the past 3 months, no extreme abnormality in meteorological conditions, and the limitation is that the influence of seismic load is not considered; the report data is reviewed and compared, the error of 3 calculation results is less than 5%, the logic audit is consistent in data transmission of each step, the expression is standard, the terms meet the geotechnical engineering investigation specification, and a complete slope stability quantitative evaluation report is generated.

[0147] In the embodiment of the present application, the strip and block division and iterative calculation ensure the safety factor calculation accuracy, the real-time pore water pressure data can dynamically reflect the influence of the rainy season on the stability, the clear threshold division can quickly determine the risk state of the sliding surface, and the problems of low efficiency and inability to adapt to dynamic environment in the traditional method are solved. The standardization processing eliminates the dimension difference of the index, the weight determined by the analytic hierarchy process conforms to the engineering practice, the clear grade division is convenient for intuitive judgment of the risk degree, and the evaluation comprehensiveness and reliability are improved.

[0148] The embodiment of the present application also provides a computing device, comprising a processor and a memory storing a computer program, wherein the computer program is run by the processor to execute the system as described above. All implementation manners in the above system embodiment are applicable to this embodiment, and the same technical effects can also be achieved.

[0149] The embodiment of the present application also provides a computer readable storage medium, which stores instructions, and when the instructions are run on a computer, the computer executes the system as described above. All implementation manners in the above system embodiment are suitable for this embodiment and can achieve the same technical effects.

[0150] The above is the preferred embodiment of the present application. It should be pointed out that for ordinary skilled in the art, without departing from the principles of the present application, several improvements and refinements can be made, which should also be considered as the protection scope of the present application.

Claims

1. A slope displacement monitoring data processing system based on UAV lidar, characterized in that, include: The acquisition module is used to acquire the raw lidar observation dataset; The processing module is used to preprocess the original lidar observation dataset to obtain an optimized lidar observation dataset. 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. 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. 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. The building module is used to perform spatial interpolation on the terrain feature dataset to construct a digital terrain surface model; 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. 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.

2. The slope displacement monitoring data processing system based on UAV lidar according to claim 1, characterized in that, The original lidar observation dataset is preprocessed to obtain an optimized lidar observation dataset, including: 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 intensity corrected by normalizing the distance and incident angle to eliminate intensity distortion caused by measurement geometry effects, generating an intensity-normalized dataset. 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 to generate a geometric correction dataset. 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.

3. The slope displacement monitoring data processing system based on UAV lidar according to claim 2, characterized in that, Based on a slope 3D observation dataset, two stable reference benchmarks are preset, an oriented 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: 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. 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; 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. 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. 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. 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.

4. The slope displacement monitoring data processing system based on UAV lidar according to claim 3, characterized in that, The parameters for land cover classification 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: 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. 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. 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. The three-dimensional point set of exposed ground surface is processed to generate an optimized terrain feature dataset.

5. The slope displacement monitoring data processing system based on UAV lidar according to claim 4, characterized in that, The terrain feature dataset is spatially interpolated to construct a digital terrain surface model, including: 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. 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. 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.

6. The slope displacement monitoring data processing system based on UAV lidar according to claim 5, characterized in that, Based on a 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: 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. 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. 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. 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.

7. The slope displacement monitoring data processing system based on UAV lidar according to claim 6, characterized in that, The displacement field calculation results are input into the risk assessment model, and the slope stability is quantitatively evaluated through a multi-source data fusion analysis platform, including: 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; 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. 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. 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.

8. The slope displacement monitoring data processing system based on UAV lidar according to claim 7, characterized in that, Based on the element 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 elements along the negative gradient direction, a set of element displacement parameters is generated, 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; 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. 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.

9. A computing device, characterized in that, include: One or more processors; A storage device for storing one or more programs, which, when executed by one or more processors, cause the one or more processors to implement the system as described in any one of claims 1 to 8.

Images