Ground surface deformation detection method based on airborne LiDAR, storage medium and equipment

Abstract

The application discloses a ground surface deformation detection method based on an airborne LiDAR, a storage medium and equipment, relates to the technical field of ground surface deformation detection, and comprises the following steps: S1, acquiring point cloud data of a detection area through an airborne LiDAR system and performing multi-period detection; S2, data preprocessing is performed to remove discrete noise points; S3, an improved progressive encryption triangulation network filtering algorithm is used to extract bare ground elevation points, and then vegetation points are further removed through vegetation index filtering, and only real ground surface points are reserved; S4, the optimized point cloud data is used to generate a digital elevation model through a Kriging interpolation method; S5, image control points of a stable area are used as reference points to test registration accuracy; S6, the digital elevation models generated in different periods are subjected to pixel-by-pixel difference calculation to obtain ground surface subsidence or uplift; and S7, a ground surface subsidence distribution map of the detection area is finally obtained. The application can realize ground surface deformation monitoring with centimeter-level precision in complex terrains and vegetation coverage areas, and has the advantages of high automation degree, high data precision and wide application range.

Description

Technical Field

[0001] This invention relates to the field of surface deformation detection technology, and in particular to a surface deformation detection method, storage medium, and device based on airborne LiDAR. Background Technology

[0002] Traditional deformation detection methods primarily employ levels, total stations, and RTK (Real-time kinematic) measurements. While these methods can achieve millimeter-level accuracy in localized areas, they suffer from drawbacks in large-scale mining area monitoring, including long data acquisition cycles, high labor costs, and difficulty in generating continuous spatial deformation fields. Furthermore, these methods only obtain linear or point-like data for cross-sections, failing to comprehensively reflect the three-dimensional variation patterns of the Earth's surface. Although synthetic aperture radar interferometry can achieve large-scale periodic coverage, it struggles to capture sudden or localized subsidence caused by coal mining, and its accuracy is limited in complex terrain areas.

[0003] With its high-density point cloud acquisition capability and centimeter-level spatial resolution, LiDAR technology is widely used in topographic mapping and surface deformation detection. Airborne LiDAR actively detects the ground surface through the time-of-flight principle of laser pulses. It can penetrate vegetation and accurately measure the elevation of bare ground, generating a high-precision digital elevation model (DEM), providing a reliable means for vertical settlement detection.

[0004] However, traditional progressive densification triangulation filtering algorithms suffer from low computational efficiency when processing large-scale point cloud data, especially in complex terrain and areas with high vegetation cover. Their filtering accuracy is also low, making it difficult to accurately extract elevation points from bare ground. Furthermore, existing semi-variogram models typically rely on fixed parameters, leading to inconsistent model accuracy under different terrain and data characteristics. Particularly for areas with significant terrain undulations or complex vegetation distribution, current technologies struggle to guarantee high-precision and high-reliability surface deformation monitoring, often neglecting the impact of local changes and failing to adapt to variable environments. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a surface deformation detection method, storage medium, and equipment based on airborne LiDAR. This method enables centimeter-level accuracy in monitoring surface deformation in complex terrain and vegetated areas. It features high automation, high data accuracy, and wide applicability, and can be widely applied in fields such as mining subsidence monitoring, landslide disaster early warning, geological environment assessment, and ecological restoration and management.

[0006] Definitions:

[0007] 1. Airborne LiDAR System: As an important payload of UAV aerial surveying platform, the airborne LiDAR system integrates multiple technologies such as laser ranging, high-precision positioning and attitude determination, and multi-source data synchronous acquisition. It mainly consists of a laser scanner, a global navigation satellite system, an inertial measurement unit, a high-resolution optical imaging module, a data storage and processing unit, and UAV platform adaptation components. Through the collaborative work of each subsystem, it achieves dynamic acquisition and real-time processing of three-dimensional spatial information of the ground surface.

[0008] 2. Digital Elevation Model (DEM): A digital simulation of ground terrain (i.e., a digital representation of the surface morphology of terrain) is achieved through limited terrain elevation data. In other words, it is a physical ground model that represents ground elevation using an ordered array of numerical values.

[0009] 3. RGB: The RGB color mode is an industry color standard that uses variations in the three color channels—red (R), green (G), and blue (B)—and their superposition to obtain a variety of colors. R represents red, G represents green, and B represents blue.

[0010] To achieve the above-mentioned technical objectives and effects, the present invention is implemented through the following technical solution:

[0011] The method for detecting land surface deformation based on airborne LiDAR includes the following steps:

[0012] S1: Acquire point cloud data of the detection area through an airborne LiDAR system and perform multi-period detection;

[0013] S2: Data preprocessing to remove discrete noise;

[0014] S3: An improved progressive densification triangular mesh filtering algorithm is used to extract bare ground elevation points, and then vegetation points are further removed by vegetation index filtering, leaving only the real surface points.

[0015] S4: Generate a digital elevation model from the optimized point cloud data using the Kriging interpolation method;

[0016] S5: Use the control points in the stable region as reference points to check the registration accuracy:

[0017] like If the accuracy is not met, return to step S1 and re-acquire the data.

[0018] like If the accuracy meets the requirements, proceed to step S6, where... The system's pre-set accuracy requirements and allowable DEM registration errors; This refers to the actual calculated DEM registration error.

[0019] S6: Perform subtraction calculations on each pixel of the digital elevation models generated at different times to obtain the amount of surface subsidence or uplift.

[0020] S7: The final result is a map showing the distribution of surface subsidence in the tested area.

[0021] Preferably, in step S2, discrete noise points are removed from the acquired point cloud data using statistical filtering:

[0022] S201, For each point Find its nearest The neighboring points are denoted as follows: ,calculate Arrive here The average Euclidean distance of the neighborhood points :

[0023] ;

[0024] S202, the average Euclidean distance from all points to their corresponding neighborhood points. Calculate the global mean and standard deviation :

[0025] ;

[0026] S203, Calculate the distance threshold :

[0027] ;

[0028] S204. Traverse all points. If the average Euclidean distance from a point to its corresponding neighborhood points is greater than the average distance from the point to the neighborhood points, then... If the value is not specified, it is marked as a noise point; otherwise, it is retained as a valid point.

[0029] in, To obtain the first Point cloud data, , The total number of points in the acquired point cloud data; It is a positive integer; is a coefficient.

[0030] Preferably, in step S3, the extraction of bare ground elevation points using the improved progressively encrypted triangular network filtering algorithm specifically includes the following steps:

[0031] S301. Determine the minimum side length of the bounding rectangle for each segment of the point cloud. and maximum points The point cloud is segmented into blocks according to the following formula:

[0032] ;

[0033] in: The side length of the bounding rectangle of the input point cloud; The total number of points in the input point cloud data; The side length of the bounding rectangle of the point cloud after segmentation;

[0034] Repeat this process until the side length of the bounding rectangle and the number of points of all segmented point clouds meet the threshold requirements.

[0035] S302. Set the sliding window size. First, set the long side of the rectangular window to be parallel to the X-axis and the short side to be parallel to the Y-axis. Starting from the origin, move the window stepwise along the X and Y axes with preset step sizes m and ε, respectively. Compare the heights of the point cloud within the window's coverage area to find the lowest point and record it as a potential ground seed point. Continue until the sliding window covers the entire point cloud area within the block. Then, adjust the long side of the rectangular window to be parallel to the Y-axis and the short side to be parallel to the X-axis, and repeat the same step-size movement rule from the origin, adding all collected points to the initial seed point set. middle;

[0036] S303. Using seed points selected three or more times as high-confidence seeds, fit a local plane to the surrounding neighborhood points with radius ζ using weighted least squares. Calculate the perpendicular distance d and slope difference Ω of each neighborhood point to the plane. If d is less than a threshold and the angle between the normal vector of the neighboring point and the normal vector of the plane is also less than a threshold, then add the point to the seed point set. middle, ;

[0037] S304, Based on the seed point set Generated using the Delaunay triangulation method , It is a triangular network. ;

[0038] in, For the seed point set The triangular mesh generated using the Delaunay triangulation method; This is the initial triangular mesh;

[0039] S305. For the remaining points that have not been added to the seed point set. Calculate the perpendicular distance from each triangle to the nearest triangle. With projection angle :

[0040] ;

[0041] ;

[0042] in, Let be the normal vector of the triangular face. Let be the projection point of the point onto the triangular face; For the first Points that have not been added to the seed point set;

[0043] If the following conditions are met:

[0044] ;

[0045] This point is then considered a ground point and added to the seed point set for the next round. ;

[0046] in, The set vertical distance threshold; The set projection angle threshold;

[0047] The seed point set is updated in each iteration:

[0048] ;

[0049] S306, If new ground points are added in this round If the iteration terminates, the output will be... That is, the set of extracted seed points; otherwise, Increment the value by 1 and re-enter step S304.

[0050] Preferably, in step S3, further removing vegetation points through vegetation index filtering specifically includes the following steps:

[0051] Calculating the visible light band vegetation index based on point cloud RGB information :

[0052] ;

[0053] ;

[0054] Set threshold ,like If so, it is identified as a vegetation point and removed;

[0055] in, Represents the red attribute Representing green attributes It represents the blue attribute.

[0056] Preferably, step S4 specifically includes the following steps:

[0057] S401, Based on the detection area range Define grid resolution And generate grid nodes: ;

[0058] in, The minimum X coordinate of the detection area; The maximum X coordinate of the detection area; The minimum Y-coordinate of the detection area; The maximum Y-coordinate of the detection area; For grid row and column indexes; The coordinates of the grid nodes to be predicted;

[0059] S402. Calculate the empirical semivariogram of ground point samples. :

[0060] ;

[0061] in, The distance between sample points; This represents the distance step size; For all sample pairs Falling into the distance range The point logarithm; For the first Elevation of each ground point; For the first Elevation of each ground point; For all points that satisfy the condition, the distance between them should fall into the range. The set of sample point pairs; It is a combination of two spatial points;

[0062] S403. Using a spherical model to fit the empirical semivariogram. :

[0063] ;

[0064] in, Value of a nugget; This is the base value; For variable range;

[0065] The fitting objective is: ;

[0066] S404, at each grid node to be predicted At this location, the elevation is estimated using ordinary kriging. :

[0067] ;

[0068] Kriging Prediction Variance The calculation is as follows:

[0069] ;

[0070] A smaller variance indicates a higher interpolation confidence level;

[0071] in, The number of ground sample points participating in the interpolation; For the first Kriging weight coefficients for each sample point; For sample points With the point to be predicted The theoretical semivariogram value corresponding to the distance between them; These are Lagrange multipliers used to introduce unbiased constraints in ordinary Kriging.

[0072] Preferably, step S4 further includes step S405:

[0073] calculate ;

[0074] in, Root mean square error is used to measure the overall error level between predicted elevation and actual elevation in cross-validation. To leave one predicted value means: in calculating the first predicted value... When predicting the elevation at the nth sample point, the nth sample point is used as the basis for the prediction. One sample point was removed from the sample set, and only the remaining sample points were used. Perform ordinary kriging interpolation on each sample point to obtain the sample points. Predicted location;

[0075] Calculate the mean error :

[0076] ;

[0077] Let the first The prediction standard deviation for each sample point is Calculate the standardized residuals:

[0078] ;

[0079] in, For the first The standardized residuals of each sample point are used to test whether the prediction error is consistent with the Kriging prediction variance.

[0080] like , The mean error threshold, and ,in, for The mean, for The variance, and If the threshold is not exceeded, the requirement is met; otherwise, the empirical semivariogram model or parameters need to be readjusted and the process proceeds to step S404 until the requirement is met.

[0081] Preferably, the formula for adjusting the parameters of the empirical semivariogram is:

[0082] ;

[0083] in, Indicates the value of the nugget Variable range or base value One of them, To adjust the direction, For the new parameters, To adjust the value size;

[0084] ;

[0085] ;

[0086] in, This is an empirical coefficient; This represents the deviation ratio. ; For target error; It is a symbolic function; The corresponding error index is as follows: when When representing the value of a nugget, This represents the difference between the actual error and the theoretical error; when When indicating a change in range, This represents the error term related to the difference between empirical error and target error; when When representing the sill value, This represents the difference between the empirical semivariogram and the model semivariogram.

[0087] Preferred, .

[0088] Correspondingly, a computer-readable storage medium stores at least one instruction, which is loaded and executed by a processor to implement the airborne LiDAR-based surface deformation detection method described in any one of the preceding claims.

[0089] Correspondingly, a computer device includes a processor and a memory, the memory storing at least one instruction, which is loaded and executed by the processor to implement the airborne LiDAR-based surface deformation detection method described in any of the above-mentioned embodiments.

[0090] The beneficial effects of this invention are:

[0091] First, this invention combines airborne laser scanning, a high-precision positioning and attitude determination system, and multi-stage data registration to extract settlement through temporal DEM differential. This not only overcomes the inefficiency and limitations of traditional measurement methods, but also enables high-precision and high-efficiency monitoring of three-dimensional surface deformation in complex terrain and vegetated areas. It can be widely applied to mining area settlement monitoring, landslide disaster early warning, and other geological environment assessment fields, providing important technical support for mining area ecological restoration, geological disaster prevention and control, and resource development.

[0092] Secondly, this invention extracts surface subsidence data using the DEM difference method, and analyzes the subsidence results and various factors affecting accuracy. Finally, it is found that the surface subsidence extracted using the method described in this invention has an accuracy of 2.5 cm, and has the advantages of high automation, high data accuracy, and wide applicability.

[0093] Third, the improved progressive densification triangulation filtering algorithm of this invention has significant advantages in extracting bare ground elevation points. The algorithm improves the processing efficiency of point cloud data through block processing and parallel computing, especially significantly reducing computation time when processing large-scale data. The improved algorithm can accurately remove vegetation interference and extract bare ground elevation points in complex terrain and high-vegetation-cover areas, solving the problem of insufficient accuracy in traditional methods. By combining vegetation index filtering, the algorithm efficiently separates ground points from non-ground points, reducing misclassification and incorrect rejection, and is particularly suitable for areas with complex vegetation. This invention not only improves efficiency and accuracy but also expands its applicability.

[0094] Fourth, the modified adjusted empirical semivariogram model or parameters have significant advantages. By introducing an adjustment mechanism, the model can dynamically optimize the semivariogram parameters according to the characteristics of the actual data, ensuring high accuracy under different terrain and data conditions. Dynamic adjustment can effectively avoid the instability of accuracy caused by fixed parameters in traditional methods. Especially in complex terrain and vegetated areas, the model can accurately reflect local terrain changes, significantly improving the accuracy and reliability of land subsidence monitoring. Attached Figure Description

[0095] Figure 1 This is a flowchart of the improved progressive encryption triangular mesh filtering algorithm for extracting bare ground elevation points according to the present invention;

[0096] Figure 2 This is a top view of the ground points obtained in an embodiment of the present invention;

[0097] Figure 3 This is a top view of a non-ground point obtained in an embodiment of the present invention;

[0098] Figure 4 This is a 45° isometric view of the ground point obtained in an embodiment of the present invention;

[0099] Figure 5 This is a 45° isometric view of a non-ground point obtained in an embodiment of the present invention;

[0100] Figure 6 This is a schematic diagram of generating a digital elevation model using the Kriging interpolation method in an embodiment of the present invention;

[0101] Figure 7 This is a schematic diagram showing the percentage of errors distributed in the range of -5mm to 5mm in an embodiment of the present invention;

[0102] Figure 8 This is a surface subsidence distribution map obtained in an embodiment of the present invention;

[0103] Figure 9 This is a schematic diagram of the statistical results of surface subsidence in an embodiment of the present invention. Detailed Implementation

[0104] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand the present invention and implement it. However, the embodiments are not intended to limit the present invention.

[0105] The method for detecting land surface deformation based on airborne LiDAR includes the following steps:

[0106] S1: Acquire point cloud data of the detection area using an airborne LiDAR system and perform multi-stage detection. This step is used to acquire data. The airborne LiDAR system is mounted on an airborne platform and performs multi-stage laser scanning of the detection area (i.e., the target area) to acquire high-density 3D point cloud data at different times. For example, high-density point cloud data can be acquired using an airborne platform equipped with a laser scanner, GPS, and IMU (Inertial Measurement Unit).

[0107] S2: Data preprocessing to remove discrete noise. This step is used to preprocess the acquired raw point cloud data. Preferably, statistical filtering is used to remove discrete noise. For example, by utilizing the average distance characteristics of each point in the point cloud to its neighboring points, combined with a Gaussian normal distribution model, outliers (i.e., isolated noise points) are detected and removed. Specifically, this includes the following steps:

[0108] S201, For each point Find its nearest The neighboring points are denoted as follows: ,calculate Arrive here The average Euclidean distance of the neighborhood points :

[0109] ;

[0110] S202, the average Euclidean distance from all points to their corresponding neighborhood points. Calculate the global mean and standard deviation :

[0111] ;

[0112] S203, Calculate the distance threshold :

[0113] ;

[0114] S204. Traverse all points. If the average Euclidean distance from a point to its corresponding neighborhood points is greater than the average distance from the point to the neighborhood points, then... If the value is not specified, it is marked as a noise point; otherwise, it is retained as a valid point.

[0115] in, To obtain the first Point cloud data, , The total number of points in the acquired point cloud data; It is a positive integer; This is a coefficient, usually 2 or 3.

[0116] S3: An improved progressive densification triangular mesh filtering algorithm is used to extract bare land elevation points, and then vegetation points are further removed by vegetation index filtering to remove vegetation interference and retain only the real surface points to obtain a high-precision point cloud of bare land area.

[0117] Preferred, such as Figure 1 As shown, the improved progressively encrypted triangular mesh filtering algorithm for extracting bare ground elevation points includes the following steps:

[0118] S301. Based on factors such as point cloud density and computer performance, determine the minimum side length of the circumscribed rectangle for each point cloud segment. and maximum points The point cloud is segmented into blocks according to the following formula:

[0119] ;

[0120] in: The side length of the bounding rectangle of the input point cloud; The total number of points in the input point cloud data; The side length of the bounding rectangle of the point cloud after segmentation;

[0121] Repeat this process until the side length and number of points of the bounding rectangle of all segmented point clouds meet the threshold requirements.

[0122] When processing large-scale point cloud data, processes such as triangulation construction, point cloud localization, and spatial relationship determination usually consume a lot of computing time. To improve processing efficiency, this invention divides the point cloud data into blocks and processes each block of point cloud data using parallel computing.

[0123] S302. Within each point cloud block, initial ground seed points are extracted using a sliding window method. First, the sliding window size is set, with the long side of the rectangular window parallel to the X-axis and the short side parallel to the Y-axis. Starting from the origin, the window moves gradually along the X and Y axes with preset step sizes m and ε, respectively, comparing the heights of the point cloud within the window's coverage area to obtain the lowest point and record it as a potential ground seed point. This process continues until the sliding window covers the entire point cloud area within the block. Subsequently, the long side of the rectangular window is adjusted to be parallel to the Y-axis and the short side to be parallel to the X-axis, and the same step size movement rule is followed to traverse the area again from the origin to extract another set of potential ground seed points. This further increases the number and spatial distribution density of initial ground seed points, adding all collected points to the initial seed point set. middle;

[0124] S303. Select seed points chosen three or more times as high-confidence seeds, and perform another round of local region growing on them. Starting from the high-confidence seed points, fit a local plane using weighted least squares method for the surrounding neighboring points with radius ζ. Calculate the perpendicular distance d and slope difference Ω of each neighboring point to the plane. If d is less than the threshold and the angle between the normal vector of the neighboring point and the normal vector of the plane is also less than the threshold, then add the point to the seed point set. middle, ;

[0125] S304, Based on the seed point set Generated using the Delaunay triangulation method , It is a triangular network. ;

[0126] in, For the seed point set The triangular mesh generated using the Delaunay triangulation method; The initial triangular mesh is generated using the starting seed point, representing a rough terrain skeleton.

[0127] S305. For the remaining points that have not been added to the seed point set. Calculate the perpendicular distance from each triangle to the nearest triangle. With projection angle :

[0128] ;

[0129] ;

[0130] in, Let be the normal vector of the triangular face. Let be the projection point of the point onto the triangular face; For the first Points that have not been added to the seed point set;

[0131] If the following conditions are met:

[0132] ;

[0133] This point is then considered a ground point and added to the seed point set for the next round. ;

[0134] in, The set vertical distance threshold; The set projection angle threshold;

[0135] The seed point set is updated in each iteration:

[0136] ;

[0137] S306, If new ground points are added in this round If the iteration terminates, the output will be... This is the extracted seed point set, which can be directly used for vegetation index filtering and DEM construction; otherwise, Increment the value by 1 and re-enter step S304.

[0138] After the point cloud data processing described above, preliminary separated ground points were obtained. However, bare ground and grass were not separated. Considering that the farmland on the mining area's surface is planted with crops, and that these crops are sown and harvested regularly, only by accurately separating the crops from the actual ground points can an accurate DEM be obtained, thus restoring the true movement of the ground surface. Preferably, further removal of vegetation points through vegetation index filtering includes the following steps:

[0139] Calculating the visible light band vegetation index based on point cloud RGB information :

[0140] ;

[0141] ;

[0142] The advantage of the visible light vegetation index is that it does not require the near-infrared band, reducing the technical requirements for sensors and enabling its widespread application to data collected by consumer drones, ordinary cameras, and other devices, allowing for the setting of threshold values. ,like If the location is identified as a vegetation point, it will be removed. Represents the red attribute Representing green attributes It represents the blue attribute.

[0143] The RGB values ​​of vegetation are affected by season and growing conditions, so it is necessary to set appropriate thresholds for each period of point cloud data. To achieve the best results, this invention sets the vegetation threshold. Set it to 0.06-0.6.

[0144] One of the main factors affecting DEM accuracy is point cloud data acquisition. Because the laser beam of an airborne LiDAR diverges with increasing distance, the density of the point cloud on the ground is uneven, especially in areas far from the flight path. This non-uniform distribution has a greater impact on planar accuracy than on elevation. As the drone's flight altitude increases, both the planar and elevation errors of the ground checkpoints also increase. However, the variation in planar error is significantly greater than that in elevation error. This is because as the drone's altitude increases, the number of points per square meter decreases, resulting in fewer points available for fitting the monitoring point's plane, thus reducing planar accuracy. For elevation, since the checkpoint is located at the center of the device's plane without undulations, the fewer points have a smaller impact on the elevation accuracy of the fitted plane's center point. Therefore, to maintain high accuracy in subsequent data, this embodiment of the invention preferentially selects a flight altitude of 78 meters as the drone's flight altitude and uses a DJI M350RTK drone equipped with an L2 LiDAR for data acquisition.

[0145] The top view of ground points obtained by the method of the present invention is as follows: Figure 2 As shown, the top view of non-ground points is as follows: Figure 3 As shown, the 45° axonometric view of the ground point is as follows. Figure 4 As shown, the 45° isometric view of non-ground points is as follows. Figure 5 As shown, where, Figure 2 and Figure 4 To the remaining ground points after removing vegetation, Figure 3 and Figure 5 As can be seen from the vegetation points to be removed, the removal effect of this invention is very good, and the separation effect is obvious. Most of the surface vegetation is removed, and it can accurately separate ground points and non-ground points with virtually no errors in removal. However, it should also be noted that voids are created after vegetation removal. Therefore, this invention needs to interpolate the voids based on the remaining real ground points.

[0146] S4: Generate a digital elevation model (DEM) from the optimized point cloud data using kriging interpolation to ensure detailed terrain representation. Specifically, step S4 includes the following steps:

[0147] S401, Based on the detection area range Define grid resolution And generate grid nodes: ;

[0148] in, The minimum X coordinate of the detection area; The maximum X coordinate of the detection area; The minimum Y-coordinate of the detection area; The maximum Y-coordinate of the detection area; For grid row and column indexes; The coordinates of the grid nodes to be predicted are shown.

[0149] Another major factor affecting DEM accuracy is the grid size during DEM generation. This invention, based on an altitude of 78m, adds control points for elevation direction optimization, generating DEMs at intervals of 0.05m, 0.1m, 0.5m, 1m, 1.5m, and 2m (i.e., grid resolution), and studies the DEM differences. Errors distributed within the range of -5mm to 5mm are statistically analyzed, and the results are as follows: Figure 7 As shown, within the range of 0.05m to 0.5m, the proportion of errors distributed within -5mm to 5mm increases with the increase of the grid size, reaching a maximum of 68.1% when the grid size is 0.5m. Further increasing the grid size reveals that the proportion of errors distributed within -5mm to 5mm begins to decrease, reaching a minimum of 59% when the grid size is 2m. Therefore, it can be concluded that when generating a DEM from point clouds, the grid resolution should be... Setting it to 0.5m will maximize the error concentration within ±5mm when performing differential processing on the resulting DEM, and will also result in a smoother differential profile curve and higher accuracy.

[0150] S402. To reflect the spatial autocorrelation of ground elevation, calculate the empirical semivariogram of ground point samples. :

[0151] ;

[0152] in, The distance between sample points (in meters); This represents the distance step size; For all sample pairs Falling into the distance range The point logarithm; For the first Elevation of each ground point; For the first Elevation of each ground point; For all points that satisfy the condition, the distance between them should fall into the range. The set of sample point pairs; This is a combination of two spatial points. To reduce computational complexity, a distance step size is typically set. m.

[0153] S403. Use a spherical model to fit the empirical semivariogram. :

[0154] ;

[0155] in, This is the nugget value (representing measurement error or microscale variation). This represents the sill value (total variance contribution). For variable range (autocorrelation distance);

[0156] The fitting objective is: Its purpose is to, given a set of empirical semivariogram values In this case, by adjusting the three parameters of the spherical model—gold value— , base value and range This makes the theoretical semivariogram calculated by the spherical model... The difference between it and the empirical semivariogram is minimal.

[0157] S404, at each grid node to be predicted At this location, the elevation is estimated using ordinary kriging. :

[0158] ;

[0159] Kriging Prediction Variance The calculation is as follows:

[0160] ;

[0161] A smaller variance indicates a higher interpolation confidence level;

[0162] in, The number of ground sample points participating in the interpolation; For the first Kriging weight coefficients for each sample point; For sample points With the point to be predicted The theoretical semivariogram value corresponding to the distance between them; These are Lagrange multipliers used to introduce unbiased constraints in ordinary Kriging.

[0163] Preferably, in DEM accuracy assessment, the root mean square error is used. To verify model stability, step S4 further includes step S405:

[0164] calculate ;

[0165] in, Root mean square error is used to measure the overall error level between predicted elevation and actual elevation in cross-validation. To leave one predicted value means: in calculating the first predicted value... When predicting the elevation at the nth sample point, the nth sample point is used as the basis for the prediction. One sample point was removed from the sample set, and only the remaining sample points were used. Perform ordinary kriging interpolation on each sample point to obtain the sample points. Predicted location.

[0166] Calculate the mean error Used to measure whether there is systematic bias in the prediction results:

[0167] ;

[0168] Let the first The prediction standard deviation for each sample point is The standardized residuals are calculated for subsequent testing.

[0169] ;

[0170] in, For the first The standardized residuals of each sample point are used to test whether the prediction error is consistent with the Kriging prediction variance.

[0171] like , The mean error threshold, ,in, ,and If the threshold is not exceeded, the requirement is met; otherwise, the empirical semivariogram model or parameters need to be readjusted and the process proceeds to step S404 until the requirement is met.

[0172] Preferably, the formula for adjusting the parameters of the empirical semivariogram is:

[0173] ;

[0174] in, Indicates the value of the nugget Variable range or base value One of them, To adjust the direction, For the new parameters (when) When representing the value of a nugget, Indicates the new nugget value; when When indicating a change in range, Indicates a new range; when When representing the sill value, (This represents the new sill value). To adjust the value size;

[0175] ;

[0176] ;

[0177] in, This is an empirical coefficient, typically ranging from 0.05 to 0.2. This represents the deviation ratio. ; For target error; It is a symbolic function; The corresponding error index is as follows: when When representing the value of a nugget, This represents the difference between the actual error and the theoretical error; when When indicating a change in range, This represents the error term related to the difference between empirical error and target error; when When representing the sill value, This represents the difference between the empirical semivariogram and the model semivariogram.

[0178] Generating a digital elevation model (DEM) using Kriging interpolation, such as... Figure 6 As shown.

[0179] S5: Use the control points in the stable region as reference points to check the registration accuracy and avoid errors caused by control point displacement.

[0180] like If the accuracy is not met, return to step S1 and re-acquire the data.

[0181] like If the accuracy meets the requirements, proceed to step S6, where... The system's pre-set accuracy requirements and allowable DEM registration errors; This represents the actual calculated DEM registration error.

[0182] S6: Perform subtraction calculations on each pixel of the digital elevation models generated at different times to obtain the amount of land subsidence or uplift, and realize the automatic extraction of the amount of land subsidence and uplift.

[0183] Since the data from different periods are all generated by Kriging interpolation, with consistent resolution and completely unified coordinate system, projection system, and datum, the DEM of Difference (DoD) model is used to calculate the change in land elevation: a change in land elevation greater than 0 indicates that the land surface is uplifted; a change in land elevation less than 0 indicates that the land surface is subsided; and a change in land elevation ≈ 0 indicates a stable zone.

[0184] S7: Finally, a surface subsidence distribution map of the detection area is obtained, which can output a surface deformation distribution map and statistical analysis results, intuitively reflecting the spatial subsidence characteristics.

[0185] In this embodiment, the data collected in the detection area were processed to obtain two data periods, DEM1 and DEM2. The surface subsidence distribution map obtained by subtracting the DEMs is shown below. Figure 8 As shown, it can be found that Figure 8 Subsidence and bulging occurred, consistent with experimental expectations, and the deformation at the test site was well reflected. Most of the ground did not experience subsidence, but errors were noted at the edges of roads and grasslands. Statistical analysis of the deformation yielded... Figure 9 ,Depend on Figure 9 It can be seen that 90.1% of the errors are distributed between -1.5cm and 1.5cm, 96.8% are distributed between -2.5cm and 2.5cm, and 63.1% are distributed between -5mm and 5mm. Therefore, the surface subsidence extraction accuracy described in this invention can be considered to be 2.5cm. This demonstrates that the DEM difference method is relatively reliable both theoretically and in terms of results.

[0186] Correspondingly, a computer-readable storage medium stores at least one instruction, which is loaded and executed by a processor to implement the airborne LiDAR-based surface deformation detection method described in any one of the preceding claims.

[0187] Correspondingly, a computer device includes a processor and a memory, the memory storing at least one instruction, which is loaded and executed by the processor to implement the airborne LiDAR-based surface deformation detection method described in any of the above-mentioned embodiments.

[0188] The above are merely preferred embodiments of the present invention and do not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made based on the content of the present invention's specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. A method for detecting surface deformation based on airborne LiDAR, characterized in that, Includes the following steps: S1: Acquire point cloud data of the detection area through an airborne LiDAR system and perform multi-period detection; S2: Data preprocessing to remove discrete noise; S3: An improved progressive densification triangular mesh filtering algorithm is used to extract bare ground elevation points, and then vegetation points are further removed by vegetation index filtering, leaving only the real surface points. S4: Generate a digital elevation model from the optimized point cloud data using the Kriging interpolation method; S5: Use the control points in the stable region as reference points to check the registration accuracy: like If the accuracy is not met, return to step S1 and re-acquire the data. like If the accuracy meets the requirements, proceed to step S6, where... The system's pre-set accuracy requirements and allowable DEM registration errors; This refers to the actual calculated DEM registration error. S6: Perform subtraction calculations on each pixel of the digital elevation models generated at different times to obtain the amount of surface subsidence or uplift. S7: The final surface subsidence distribution map of the tested area is obtained; Step S3, which uses an improved progressively encrypted triangular mesh filtering algorithm to extract bare ground elevation points, specifically includes the following steps: S301. Determine the minimum side length of the bounding rectangle for each segment of the point cloud. and maximum points The point cloud is segmented into blocks according to the following formula: ； in: The side length of the bounding rectangle of the input point cloud; The total number of points in the input point cloud data; The side length of the bounding rectangle of the point cloud after segmentation; Repeat this process until the side length of the bounding rectangle and the number of points of all segmented point clouds meet the threshold requirements. S302. Set the sliding window size. First, set the long side of the rectangular window to be parallel to the X-axis and the short side to be parallel to the Y-axis. Starting from the origin, move the window stepwise along the X and Y axes with preset step sizes m and ε, respectively. Compare the heights of the point cloud within the window's coverage area to find the lowest point and record it as a potential ground seed point. Continue until the sliding window covers the entire point cloud area within the block. Then, adjust the long side of the rectangular window to be parallel to the Y-axis and the short side to be parallel to the X-axis, and repeat the same step size movement rule starting from the origin, adding all collected points to the initial seed point set. middle; S303. Using the data collection points obtained three or more times as high-confidence seeds, fit a local plane to the surrounding neighborhood points with a radius of ζ using the weighted least squares method. Calculate the perpendicular distance d and slope difference Ω of each neighborhood point to the plane. If d is less than a threshold and the angle between the normal vector of the neighboring point and the normal vector of the plane is also less than a threshold, then add the point to the seed point set. middle, ; S304, Based on the seed point set Generated using the Delaunay triangulation method , It is a triangular network. ; in, For the seed point set The triangular mesh generated using the Delaunay triangulation method; This is the initial triangular mesh; S305. For the remaining points that have not been added to the seed point set. Calculate the perpendicular distance from each triangle to the nearest triangle. With projection angle : ； ； in, Let be the normal vector of the triangular face. Let be the projection point of the point onto the triangular face; For the first Points that have not been added to the seed point set; If the following conditions are met: ； This point is then considered a ground point and added to the seed point set for the next round. ; in, The set vertical distance threshold; The set projection angle threshold; The seed point set is updated in each iteration: ； S306, If new ground points are added in this round If the iteration terminates, the output will be... That is, the set of extracted seed points; otherwise, Increment the value by 1 and re-enter step S304.

2. The surface deformation detection method based on airborne LiDAR according to claim 1, characterized in that, In step S2, the acquired point cloud data is processed using statistical filtering to remove discrete noise points. S201, For each point Find its nearest The neighboring points are denoted as follows: ,calculate Arrive here The average Euclidean distance of the neighborhood points : ； S202, the average Euclidean distance from all points to their corresponding neighborhood points. Calculate the global mean and standard deviation : ； S203, Calculate the distance threshold : ； S204. Traverse all points. If the average Euclidean distance from a point to its corresponding neighborhood points is greater than the average distance from the point to the neighborhood points, then... If the value is not specified, it is marked as a noise point; otherwise, it is retained as a valid point. in, To obtain the first Point cloud data, , The total number of points in the acquired point cloud data; It is a positive integer; is a coefficient.

3. The surface deformation detection method based on airborne LiDAR according to claim 1, characterized in that, In step S3, the further removal of vegetation points through vegetation index filtering specifically includes the following steps: Calculating the visible light band vegetation index based on point cloud RGB information : ； ； Set threshold ,like If so, it is identified as a vegetation point and removed; in, Represents the red attribute Representing green attributes It represents the blue attribute.

4. The surface deformation detection method based on airborne LiDAR according to claim 3, characterized in that, Step S4 specifically includes the following steps: S401, Based on the detection area range Define grid resolution And generate grid nodes: ; in, The minimum X coordinate of the detection area; The maximum X coordinate of the detection area; The minimum Y-coordinate of the detection area; The maximum Y-coordinate of the detection area; For grid row and column indexes; The coordinates of the grid nodes to be predicted; S402. Calculate the empirical semivariogram of ground point samples. : ； in, The distance between sample points; This represents the distance step size; For all sample pairs Falling into the distance range The point logarithm; For the first Elevation of each ground point; For the first Elevation of each ground point; For all points that satisfy the condition, the distance between them should fall into the range. The set of sample point pairs; It is a combination of two spatial points; S403. Using a spherical model to fit the empirical semivariogram. : ； in, Value of a nugget; This is the base value; For variable range; The fitting objective is: ; S404, at each grid node to be predicted At this location, the elevation is estimated using ordinary kriging. : ； Kriging Prediction Variance The calculation is as follows: ； A smaller variance indicates a higher interpolation confidence level; in, The number of ground sample points participating in the interpolation; For the first Kriging weight coefficients for each sample point; For sample points With the point to be predicted The theoretical semivariogram value corresponding to the distance between them; These are Lagrange multipliers used to introduce unbiased constraints in ordinary Kriging.

5. The surface deformation detection method based on airborne LiDAR according to claim 4, characterized in that, Step S4 also includes step S405: calculate ; in, Root mean square error is used to measure the overall error level between predicted elevation and actual elevation in cross-validation. To leave one predicted value means: in calculating the first predicted value... When predicting the elevation at the nth sample point, the nth sample point is used as the basis for the prediction. One sample point was removed from the sample set, and only the remaining sample points were used. Perform ordinary kriging interpolation on each sample point to obtain the sample points. Predicted location; Calculate the mean error : ； Let the first The prediction standard deviation for each sample point is Calculate the standardized residuals: ； in, For the first The standardized residuals of each sample point are used to test whether the prediction error is consistent with the Kriging prediction variance. like , The mean error threshold, and ,in, for The mean, for The variance, and If the threshold is not exceeded, the requirement is met; otherwise, the empirical semivariogram model or parameters need to be readjusted before proceeding to step S404 until the requirement is met.

6. The surface deformation detection method based on airborne LiDAR according to claim 5, characterized in that, The formula for adjusting the parameters of the empirical semivariogram is: ； in, Indicates the value of the nugget Variable range or base value One of them, To adjust the direction, For the new parameters, To adjust the value size; ； ； in, This is an empirical coefficient; This represents the deviation ratio. ; For target error; It is a symbolic function; The corresponding error index is as follows: when When representing the value of a nugget, This represents the difference between the actual error and the theoretical error; when When indicating a change in range, This represents the error term related to the difference between empirical error and target error; when When representing the sill value, This represents the difference between the empirical semivariogram and the model semivariogram.

7. The surface deformation detection method based on airborne LiDAR according to claim 4, characterized in that, 。 8. A computer-readable storage medium, characterized in that, The storage medium stores at least one instruction, which is loaded and executed by a processor to implement the airborne LiDAR-based surface deformation detection method as described in any one of claims 1-7.

9. A computer device, characterized in that, The computer device includes a processor and a memory, the memory storing at least one instruction, which is loaded and executed by the processor to implement the airborne LiDAR-based surface deformation detection method as described in any one of claims 1-7.

Images