Surface subsidence monitoring method and system based on laser point cloud data

Abstract

The present invention discloses a surface subsidence monitoring method and system based on laser point cloud data. The method comprises acquiring first-phase and second-phase point cloud data of a target area at different time periods; determining the Z-coordinate deviation of each flight path in each phase of the point cloud data; determining the Z-coordinate correction for each flight path based on the flight path deviation, and correcting the corresponding flight path using the Z-coordinate correction; and determining the surface subsidence of the target area based on the corrected flight paths using the first-phase and second-phase point cloud data. The present invention can improve the relative accuracy between point cloud data flight paths, enabling differential calculation of subsidence, resulting in more accurate subsidence monitoring results.

Description

Technical Field

[0001] The invention discloses a surface subsidence monitoring method and system based on laser point cloud data, belonging to the technical field of subsidence monitoring. Background Art

[0002] Surface subsidence monitoring is fundamental to geological disaster prevention and early warning. Traditional surface subsidence monitoring is limited by its small scope, low point density, and long monitoring cycles, limiting its ability to monitor the macroscopic characteristics and spatiotemporal evolution of geological hazards.

[0003] Airborne LiDAR (LiDAR) can rapidly acquire large-area laser point cloud data, providing a true three-dimensional representation of the Earth's surface. Therefore, it can be used for surface subsidence monitoring. However, the three-dimensional coordinates of LiDAR-derived point clouds often contain errors. These errors arise from a variety of sources, such as laser ranging errors, sensor installation errors, and position and positioning system (POS) errors. Consequently, the surface subsidence determined using LiDAR point cloud data is inaccurate.

[0004] To reduce the three-dimensional coordinate errors of point clouds, existing technologies typically use LiDAR Strip Adjustment (LSA) to eliminate errors. However, most state-of-the-art LSA methods assume sufficient POS accuracy. However, this is only an ideal state. Due to the wide range of airborne LiDAR, the overall LSA method has difficulty focusing on a specific study area and does not account for systematic errors. This results in large errors in point cloud coordinates, and its accuracy cannot meet the requirements of surface subsidence monitoring. Summary of the Invention

[0005] The purpose of this application is to provide a surface settlement monitoring method and system based on laser point cloud data, so as to solve the technical problem of inaccurate monitoring results due to large errors in point cloud coordinates in existing surface deformation monitoring.

[0006] A first aspect of the present invention provides a surface subsidence monitoring method based on laser point cloud data, comprising:

[0007] Step 1: Obtain the first phase point cloud data and the second phase point cloud data of the area to be measured in different time periods;

[0008] Step 2: Determine the flight path deviation of each flight path in the Z coordinate direction in each phase of point cloud data;

[0009] Step 3: Determine the Z coordinate correction value of each flight strip according to the flight strip deviation, and use the Z coordinate correction value to correct the corresponding flight strip;

[0010] Step 4: Determine the surface settlement of the area to be measured based on the first phase point cloud data and the second phase point cloud data after the corrected flight path.

[0011] Preferably, the step 2 specifically includes:

[0012] Step 21: Obtain the overlapping portion of two adjacent flight strips in each phase of point cloud data, and record each point in the upper flight strip in the overlapping portion as P i , each point in the lower flight path is denoted as P j ;

[0013] Step 22: Determine each P i and the nearest multiple P j The distance between the fitting planes is recorded as P i The point cloud deviation value of the point;

[0014] Step 23: Based on multiple P i The point cloud deviation value of the point and the corresponding X-axis coordinate value are used to determine the flight path deviation.

[0015] Preferably, the step 23 specifically includes:

[0016] According to multiple P i The point cloud deviation value and the corresponding X-axis coordinate value of the point are used to determine the flight path deviation constant using the polynomial fitting method;

[0017] The flight path deviation of each flight path in the Z coordinate direction in each period of point cloud data is determined according to the flight path deviation constant.

[0018] Preferably, determining the Z coordinate correction value of each flight strip according to the flight strip deviation specifically includes:

[0019] The Z coordinate correction value of each flight strip is determined according to the first formula, which is:

[0020]

[0021] Where V is the deviation error of each flight strip, A is the coefficient matrix, is the Z coordinate correction number of each flight strip, and l is the expected value of the flight strip deviation.

[0022] Preferably, the Z coordinate correction number Including Z coordinate correction of stable flight zone and Z coordinate correction of unstable flight zone but and Determined according to the second formula, the second formula is:

[0023]

[0024] Where M = N 22 -N 21 N 11- N 12 , is the transpose of A1, - is the sign of the inverse moment, P is the weight matrix, l is the expected value of the flight path deviation, and A1 and A2 are the matrices after the coefficient matrix A is decomposed.

[0025] Preferably, the step 4 specifically includes:

[0026] The surface settlement of the area to be measured is determined using the M3C2 algorithm based on the first phase point cloud data and the second phase point cloud data after the corrected flight path.

[0027] Preferably, after step 1, the method further comprises:

[0028] Eliminating outliers in the first-phase point cloud data and the second-phase point cloud data;

[0029] Resampling the first phase point cloud data and the second phase point cloud data after removing outliers using an octree method;

[0030] Accordingly, step 2 is as follows:

[0031] Determine the flight path deviation of each flight path in the Z coordinate direction in each phase of the resampled point cloud data.

[0032] Preferably, after resampling the first phase point cloud data and the second phase point cloud data after removing outliers by using an octree method, the method further includes:

[0033] Obtain the point source ID of the midpoint of the first phase point cloud data and the second phase point cloud data after resampling;

[0034] The first phase point cloud data and the second phase point cloud data are respectively split into flight strips according to the point source ID.

[0035] Preferably, after the first phase point cloud data and the second phase point cloud data are respectively split into flight strips according to the point source ID, the method further includes:

[0036] Using a point cloud ground point filtering method, the first phase point cloud data and the second phase point cloud data after the flight strip splitting are filtered respectively to obtain ground points in the first phase point cloud data and ground points in the second phase point cloud data;

[0037] Accordingly, the step 2 is specifically as follows:

[0038] Determine the flight path deviation of each flight path in the Z coordinate direction at each ground point.

[0039] A second aspect of the present invention provides a surface subsidence monitoring system based on laser point cloud data, comprising:

[0040] A data acquisition module, the data acquisition module is used to acquire the first phase point cloud data and the second phase point cloud data of the area to be measured in different time periods;

[0041] A deviation determination module, the deviation determination module is used to determine the flight path deviation of each flight path in each period of point cloud data in the Z coordinate direction;

[0042] a flight strip correction module, the flight strip correction module being used to determine a Z coordinate correction number for each flight strip according to the flight strip deviation, and to correct the corresponding flight strip using the Z coordinate correction number;

[0043] A settlement determination module is used to determine the surface settlement of the area to be measured based on the first phase point cloud data and the second phase point cloud data after the corrected flight path.

[0044] Compared with the existing technology, the surface subsidence monitoring method and system based on laser point cloud data of the present invention have the following beneficial effects:

[0045] The method and system of the present invention can improve the relative accuracy of point cloud data between flight strips without flight quality data, so that differential calculation of settlement can be performed, and the surface settlement monitoring results obtained are more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 Schematic diagram of the process of a surface subsidence monitoring method based on laser point cloud data in an embodiment of the present invention;

[0047] Figure 2 Schematic diagram of the structure of a surface subsidence monitoring system based on laser point cloud data in an embodiment of the present invention;

[0048] Figure 3 Comparison of surface difference results using different methods, where (a) is the unprocessed point cloud difference result, (b) is the point cloud difference result of adjacent flight strip registration using the ICP method, and (c) is the result obtained using the quasi-stabilization adjustment method of the present invention.

[0049] Figure 4 To use the section line a Figure 3 Result diagrams of profile analysis, where (a) is the profile diagram corresponding to the unprocessed point cloud difference result, (b) is the profile diagram corresponding to the point cloud difference result of adjacent flight strip registration using the ICP method, and (c) is the profile diagram corresponding to the result obtained by the quasi-stabilization adjustment of the present invention;

[0050] Figure 5 To use the section line b Figure 3Result diagrams of profile analysis, where (a) is the profile diagram corresponding to the unprocessed point cloud difference result, (b) is the profile diagram corresponding to the point cloud difference result of adjacent flight strip registration using the ICP method, and (c) is the profile diagram corresponding to the result obtained by the quasi-stabilization adjustment of the present invention;

[0051] Figure 6 for Figure 4 Comparison of Savitzky-Golay experimental results along the middle section line a;

[0052] Figure 7 for Figure 5 Comparison of Savitzky-Golay experimental results along the middle section line b;

[0053] Figure 8 It is the settlement profile position map;

[0054] Figure 9 This is a comparison chart of sedimentation rates;

[0055] Figure 10 This is a comparison chart of the flight zone deviation before and after the adjustment in 2009;

[0056] Figure 11 This is a comparison chart of the flight zone deviation before and after the adjustment in 2016.

[0057] In the figure, 101 is a data acquisition module; 102 is a deviation determination module; 103 is a flight strip correction module; and 104 is a settlement determination module. DETAILED DESCRIPTION

[0058] In the following description, specific details such as particular system structures and techniques are provided for purposes of illustration, not limitation, to facilitate a thorough understanding of the embodiments of the present invention. However, it will be apparent to those skilled in the art that the present invention may be practiced in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted so as not to obscure the description of the present invention with unnecessary detail.

[0059] The first aspect of the present invention provides a surface subsidence monitoring method based on laser point cloud data, such as Figure 1 Shown, including:

[0060] Step 1: Obtain the first phase point cloud data and the second phase point cloud data of the area to be measured in different time periods.

[0061] In the embodiment of the present invention, both the first phase point cloud data and the second phase point cloud data are original point cloud data.

[0062] Prior art methods for improving point cloud data accuracy can generally be categorized as sensor-driven and data-driven. Sensor-driven methods involve correcting point cloud data by recording aircraft runtime status information and sensor connection methods, combined with various parameter information. This method requires high-quality flight parameter information, which is sometimes difficult to obtain. Therefore, to improve the universality of point cloud processing methods, the present invention utilizes a data-driven approach, directly starting from the point cloud data itself and combining its information record files with relevant algorithms to improve data accuracy. The specific process is as follows:

[0063] Step a: remove outliers from the first phase point cloud data and the second phase point cloud data.

[0064] Laser scanning typically produces point cloud datasets with uneven density. Furthermore, sparse outliers are generated during the measurement process. To ensure the accuracy of subsequent calculations of track deviations, this paper removes outliers from the point cloud data using the following method:

[0065] When analyzing each point, a preset number of adjacent points are selected. The number of adjacent points can be set as needed, for example, 30, 40, 50, etc. In the embodiment of the present invention, 40 adjacent points are selected and the standard deviation multiple is set to 1. When the distance of a point exceeds the average distance by more than one standard deviation, the point is marked as an outlier and will be removed.

[0066] Because the flight paths in the point cloud data have overlapping parts, resulting in uneven point cloud density, the embodiment of the present invention further provides step b to obtain a point cloud with uniform density.

[0067] Step b: resample the first phase point cloud data and the second phase point cloud data after removing outliers using the octree method, specifically:

[0068] The point cloud is resampled by replacing all points in each cell of the octree (at a certain subdivision level) with its centroid. This can control the number of output point cloud points and obtain the first and second phase point cloud data with uniform density after resampling.

[0069] Many open data sources in the prior art provide a large amount of laser point cloud data. These data are all spliced ​​whole point cloud data. Since the first phase point cloud data and the second phase point cloud data in the embodiment of the present invention are both raw data, they need to be split, specifically as follows:

[0070] Step c: Obtain the point source IDs of the points in the first and second point cloud data phases. Then, split the first and second point cloud data phases into flight zones based on the point source IDs. Taking the laz point cloud format as an example, the industry standard format for LiDAR data is a binary file, as shown in Table 1. This embodiment of the present invention can split the point cloud into flight zones based on the point source ID information, and store the point cloud data for each flight zone separately.

[0071] Table 1 Point cloud data information record table

[0072] Item Format Size X Long 4 bytes Y Long 4 bytes Z Long 4 bytes Intensity Unsigned short 2 bytes Return Number 3 bits (bits0,1,2) 3 bits Scan Direction 3 bits (bits3,4,5) 3 bits Edge of Flight Line 1 bit (bit6) 1 bit Classification 1 bit (bit7) 1 bit Scan Ange Rank Unsigned char 1 byte User Data Unsigned char 1 byte Point Source ID Unsigned char 2 bytes

[0073] To further improve the accuracy of the data, the embodiment of the present invention also filters the first phase point cloud data and the second phase point cloud data, specifically:

[0074] Step d: Filter the first phase point cloud data and the second phase point cloud data respectively using the point cloud ground point filtering method to obtain the ground points in the first phase point cloud data and the ground points in the second phase point cloud data.

[0075] The point cloud ground point filtering method in the embodiment of the present invention is an existing filtering method, which is also called cloth filtering (cerebrospinal method). This method is based on cloth simulation and is a 3D computer graphics algorithm used to simulate cloth in a computer program. First, a basic formula (1) is defined:

[0076]

[0077] Where m is the weight of the particle, X represents the position of the particle in the "cloth" at time t, and F ext (X, t) represents external driving factors (gravity, collision, etc.), F int (X,t) represents the internal driving factor (internal connection between particles). Assume that there is only external factor F ext (X,t), and set the internal factor F int (X, t) is 0, and we can get formula (2):

[0078]

[0079] Where m is the weight of the particle, usually set to 1, and △t is the time step. Since G is a constant, the new position of the particle in the next iteration can be calculated as long as △t is given. In order to constrain the inversion problem of the particle in the blank area of ​​the inversion surface, the internal factor F is considered. int(X, t). Randomly select two adjacent particles. If both particles are movable, move them the same distance in opposite directions. If one particle is immovable, move the other particle. If both particles have the same height, do not move them. The displacement can be calculated using formula (3):

[0080]

[0081] Where d is the displacement of the particle; when the particle is movable, b = 1, and when it is immovable, b = 0; p i is the neighboring particle of p0, and n is the unit vector (0, 0, 1) that normalizes the point to the vertical direction. Therefore, the brief steps of cloth simulation are:

[0082] 1) Invert the denoised and resampled lidar point cloud.

[0083] 2) Initial cloth grid, set the grid size (gridresolution, gr). The initial "cloth" position is usually above the highest point.

[0084] 4) Project all lidar points and grid particles onto the same horizontal plane, find the nearest neighbor point (corresponding point, cp) of each particle, and record its elevation before projection (intersection height value, ihv)

[0085] 5) For each movable grid "particle", calculate its displacement caused by gravity and compare it with the elevation of the current particle's corresponding nearby point. If the particle's height is lower than or equal to the pre-projection elevation value, set the particle's height to ihv and set it as an immovable point.

[0086] 6) For each grid “particle”, calculate the displacement caused by the internal driving factors.

[0087] 7) Repeat steps 5) and 6) above until the maximum elevation change of all particles is small enough or the number of iterations reaches a preset value, then stop the simulation process.

[0088] 8) Calculate the height difference between the lidar point cloud and the grid particles to distinguish between ground points and non-ground points. If the distance between the laser point and the simulated particle is less than a preset threshold, it is considered to be a ground point, otherwise it is considered to be a non-ground point.

[0089] The order of step c and step d in the embodiment of the present invention can be changed.

[0090] Since there are systematic deviations in the z-direction coordinate values ​​between the flight strips, the embodiment of the present invention adopts a flight strip quasi-adjustment method to eliminate the errors, as shown in the following steps 2 and 3.

[0091] Step 2: Determine the flight path deviation of each flight path in the Z coordinate direction in each period of point cloud data.

[0092] After obtaining the ground point, the embodiment of the present invention calculates the flight path deviation as follows:

[0093] Determine the flight path deviation of each flight path in the Z coordinate direction at each ground point.

[0094] The above step 2 specifically includes:

[0095] Step 21: Obtain the overlapping portion of two adjacent flight strips in each phase of point cloud data, and record each point in the upper flight strip in the overlapping portion as P i , each point in the lower flight path is denoted as P j .

[0096] In one embodiment, the point cloud data used in this step is point cloud data of ground points.

[0097] Step 22: Determine each P i and the nearest multiple P j The distance between the fitting planes is recorded as P i The point cloud deviation value of the point.

[0098] For example, using P i The 6 closest points j Point fitting plane, calculate P i The distance from the point to the fitting plane is taken as P i The point cloud deviation value D of the point.

[0099] Step 23: Based on multiple P i The point cloud deviation value and the corresponding x-axis coordinate value of the point are used to determine the flight path deviation, including:

[0100] Step 231: Based on multiple P i The point cloud deviation value and the corresponding x-axis coordinate value of the point are used to determine the flight path deviation constant z0 using the polynomial fitting method, where the fitting formula is as follows:

[0101] D=z0+bx (4)

[0102] b in formula (4) is the fitting constant.

[0103] Step 232: Determine the flight path deviation of the corresponding flight path according to the flight path deviation constant z0.

[0104] In an embodiment of the present invention, steps 231 and 232 are repeated to obtain multiple deviation constants z0 for the same flight strip, and their average is taken, and the average is used as the flight strip deviation of the corresponding flight strip. In an embodiment of the present invention, the number of repetitions can be 5 times, 10 times, or 20 times, etc.

[0105] Step 3: According to the flight strip deviation, use the adjustment model to determine the Z coordinate correction number of each flight strip, and use the Z coordinate correction number to correct the corresponding flight strip.

[0106] In the embodiment of the present invention, the Z coordinate correction number of each flight strip is determined according to formula (5):

[0107]

[0108] Where V is the deviation error of each flight strip, A is the coefficient matrix, is the Z coordinate correction number of each flight strip, l is the expected value of the flight strip deviation, which is equal to the mean of multiple deviation constants z0.

[0109] For example, the first phase point cloud data and the second phase point cloud data in the present invention are split into 14 flight strips. Through the comparative analysis of GNSS and surface continuity, the 7th and 8th flight strips in the 14 flight strips are determined to be stable flight strips. Therefore, the coefficient matrix A can be decomposed into A1 and A2. Let represents the correction number of the unknown parameters of the stable point, represents the correction number of the unstable point parameter, and the error can be written as:

[0110]

[0111] in:

[0112]

[0113]

[0114]

[0115]

[0116] Construct the normal equation:

[0117]

[0118] in,

[0119] Eliminate unstable unknown parameter corrections using reduction methods The correction number of unknown parameters containing only stable points can be obtained now that:

[0120]

[0121] in,

[0122] Because there is no starting data, the minimum norm condition of the stable flight zone must be added as a constraint:

[0123]

[0124] According to the generalized inverse theory, the minimum norm solution is obtained as:

[0125]

[0126] Where, is the Z coordinate correction number for the stable flight zone, is the generalized inverse matrix of M.

[0127] Substituting it into the normal equation, we get:

[0128]

[0129] is the Z coordinate correction number of the unstable flight zone, - is the inverse matrix operation symbol, is the transpose of A1, P is the measurement weight matrix, l is the expected value of the flight path deviation, and its value is equal to the average value of multiple deviations of the same flight path. A1 and A2 are matrices after the coefficient matrix A is decomposed.

[0130] After obtaining the Z coordinate correction number of each flight strip, use the Z coordinate correction number to correct the corresponding flight strip, specifically:

[0131] The corrected point cloud is obtained by adding the Z coordinate of each flight path point cloud to its respective correction number.

[0132] Then perform accuracy assessment:

[0133] According to formula (11), the unit weighted mean error σ0 of the correction number can be obtained.

[0134]

[0135] Where P is the measurement weight matrix, V is the correction number, and r is the number of redundant observations.

[0136] From the cofactor propagation rate, the cofactor matrix of unknown parameters can be obtained as Q ZZ , using formula (12), we can get the mean error σ of the flight strip correction number: zi for:

[0137]

[0138] The co-factor of the flight path can be obtained from the co-factor propagation rate: φφ=FQ zz F T , and finally according to the formula The mean error σ of the flight path can be obtained φ0 .

[0139] Step 4: Based on the first-phase point cloud data and the second-phase point cloud data after the corrected flight path, the surface settlement of the area to be measured is determined using the M3C2 algorithm.

[0140] To obtain the deformation information of the two point clouds, the M3C2 point cloud comparison algorithm, first proposed by Lague, is used. This algorithm can detect changes in complex terrain directly on the point cloud without the need for meshing, and the change calculation is less affected by spatial point density, surface roughness, and different sampling positions. The main steps are:

[0141] 1) Select the core point cloud. Using the core point cloud instead of the entire point cloud as the main body of the calculation can reduce complexity, save time and improve calculation efficiency.

[0142] 2) Calculate the 3D surface normal. Set the normal scale to determine the diameter of the spherical neighborhood around each core point used to calculate the local normal. This normal is used to locate the cylinder within which the equivalent point in the other cloud is searched.

[0143] 3) Point cloud distance calculation. All core point clouds are traversed to obtain shape information. Confidence intervals for spatial variables are determined. The M3C2 algorithm uses the positional variability along N within these points (i.e., the local roughness in the normal direction) as a measure of the uncertainty of their average position, thereby determining the confidence interval for the distance measurement.

[0144] The method of this application starts from the perspective of the data itself, and improves the accuracy of point cloud data through flight strip adjustment, so that it can perform differential calculation of deformation variables.

[0145] The second aspect of the present invention provides a surface subsidence monitoring system based on laser point cloud data, such as Figure 2 As shown, it includes a data acquisition module 101 , a deviation determination module 102 , a flight strip correction module 103 and a settlement determination module 104 .

[0146] The data acquisition module 101 is used to acquire the first phase point cloud data and the second phase point cloud data of the area to be measured in different time periods;

[0147] The deviation determination module 102 is used to determine the flight path deviation of each flight path in the Z coordinate direction in each period of point cloud data;

[0148] The flight strip correction module 103 is used to determine the Z coordinate correction value of each flight strip according to the flight strip deviation, and use the Z coordinate correction value to correct the corresponding flight strip;

[0149] The settlement determination module 104 is used to determine the surface settlement of the area to be measured based on the first phase point cloud data and the second phase point cloud data after the corrected flight path.

[0150] The method of the present invention will be described in detail below with more specific examples.

[0151] 1. Comparative Analysis

[0152] 1.1 Surface area comparison analysis

[0153] The point cloud data of the same area in 2009 and 2016 were obtained from an open source database. The point cloud data in the open source database has been subjected to ALS point cloud adjustment for a large area. However, in order to take into account the absolute accuracy of the overall point cloud, the post-processing accuracy of the point cloud flight strips after adjustment is 10cm. When focusing on a specific study area, the deviation between flight strips cannot meet the needs of settlement monitoring. In order to obtain a more continuous and accurate deformation field, the point cloud of the study area needs to be refined. The embodiment of the present invention adopts different methods to obtain the deformation results of the study area. Figure 3 As shown, Figure 3 (a) to (c) in the figure are respectively the unprocessed point cloud difference results, the point cloud difference results of adjacent flight strip registration using the Iterative Closest Point method (ICP method for short), and the point cloud difference results obtained by the Quasi-Stable Adjustment (QSA) of the present invention. Figure 3 It can be seen that the deformation of the unprocessed point cloud perpendicular to the heading direction shows discontinuous changes, which is caused by the large error of the flight strip. From the point cloud difference results of the adjacent flight strips registered by the ICP method, it can be seen that the deformation of the middle deformation area perpendicular to the heading direction is relatively continuous, and the deformation changes at both ends fluctuate greatly. The ICP algorithm only aligns the point clouds between adjacent flight strips. The final registered flight strip will be affected by the cumulative error, which reduces the absolute accuracy of the point cloud. Since the premise of differential calculation of settlement is to have two periods of absolutely high point cloud data, the ICP method causes the coordinates of the final registered flight strip point cloud to change greatly. Therefore, the difference results of the two-phase point cloud will fluctuate greatly on both sides. From the domain comparison diagram, it can be seen that the differential results obtained by the quasi-stable adjustment method proposed by the present invention are relatively continuous in the deformation changes in the heading direction and the direction perpendicular to it.

[0154] 1.2 Comparative analysis of section lines

[0155] In order to better verify the effectiveness of the method of the present invention, the above domain results are analyzed as follows Figure 4 and Figure 5 shown.

[0156] Compared with other data, point cloud data is relatively disordered and has chaotic elevation jumps. Since the line graph is greatly affected by a single jump point, the Savitzky-Golay algorithm is used to remove elevation mutation noise in the result to better reflect the trend information. Figure 6 、 Figure 7 shown.

[0157] The number of extreme points in the settlement curve can reflect the quality of settlement monitoring. The more extreme points there are, the worse the continuity of deformation on the surface domain. The dotted line represents the experimental results of the unprocessed raw data, which has a large number of extreme points. The black solid curve and the gray solid curve represent the experimental results of the ICP algorithm and the method of the present invention, respectively. It is not difficult to see that the results are highly consistent in the middle part, and the ICP algorithm has large fluctuations at both ends of the curve. This is due to the limitations of the ICP algorithm. When adjacent flight strips are aligned, the cumulative error increases with the number of alignments, which is more obvious in the surface domain.

[0158] 2. Single point verification

[0159] Use Figure 8 The deformation rates of the two GNSS points in the figure verify the accuracy of the settlement map obtained by the method of the present invention. The settlement rate is as follows: Figure 9 shown.

[0160] The sedimentation rate at the LRCA GNSS point is -2.75 cm / year, and the sedimentation rate monitored by lidar is -2.52 cm / year. The sedimentation rate at the ORCA GNSS point is -8.42 cm / year, and the sedimentation rate monitored by lidar is -7.25 cm / year. The results are relatively close.

[0161] 3. Accuracy assessment

[0162] After the flight strip quasi-adjustment calculation, the flight strip correction numbers are shown in Table 2:

[0163] Table 2 Z coordinate corrections for flight strips

[0164]

[0165] According to the formula The unit weighted mean error σ0 of the corrections in 2009 and 2016 is 0.4462cm and 0.4617cm respectively. The co-factor matrix of the unknown parameters can be obtained from the co-factor propagation rate as Q ZZ09 and Q ZZ10 , using the formula The error in the obtained flight strip correction number is shown in Table 3:

[0166] Table 3 Mean error of flight strip corrections

[0167]

[0168] By the cofactor propagation rate σ φφ =FQ zz F T , The mean error of the track in 2009 and 2016 is 1.208cm and 1.190cm.

[0169] After correction, the track deviations were re-measured and shown in Tables 4 and 5.

[0170] Table 4 Statistics of track deviations after 2009 quasi-stabilization correction

[0171] flight path 1-2 2-3 3-4 4-5 5-6 6-7 7-8 8-9 9-10 10-11 11-12 12-13 13-14 <![CDATA[RMSE track ]]> 1 1.651 -1.0298 3.3448 1.0702 -1.5867 -2.7481 1.0601 -3.82 1.5547 2.1345 2.471 -3.3295 -1.1776 2 1.2793 -1.1766 3.1956 1.4888 -1.8027 -2.6936 0.8734 -3.7729 1.695 2.0501 2.328 -2.9505 -1.2029 3 1.0332 -0.7602 3.2739 1.1949 -1.6933 -2.5903 0.9585 -3.8068 1.7955 1.909 2.226 -3.4481 -0.8808 4 1.2397 -0.8424 3.4283 1.1477 -1.5907 -2.7341 0.9783 -3.6122 1.7297 2.0858 2.4198 -3.2243 -0.6693 5 1.5152 -0.7705 3.5306 1.2056 -1.5569 -2.899 0.7782 -3.6329 1.5252 1.8611 2.2663 -3.4976 -0.6424 6 1.5563 -1.129 3.4325 1.2254 -1.5344 -2.7731 0.981 -3.7849 1.6143 2.1356 2.2769 -3.3455 -1.1452 7 1.3814 -0.8961 3.0112 1.3279 -1.9046 -2.786 1.0474 -3.6531 1.9171 2.0998 2.34 -3.4991 -0.9621 8 1.4804 -0.9014 3.4049 1.3372 -1.8707 -2.517 0.9132 -3.5274 1.5078 1.828 2.1836 -3.3676 -1.0368 9 1.6255 -1.2263 3.2631 1.2375 -1.5511 -2.69 1.1181 -3.6079 1.9319 1.8282 2.1196 -3.4003 -0.9334 10 1.0384 -1.0556 3.2627 1.2389 -1.7033 -2.6523 1.0657 -3.6472 1.539 1.8763 2.2895 -3.3123 -1.1146 <![CDATA[ vz ]]> 1.38 -0.9788 3.3148 1.2474 -1.6794 -2.7083 0.9774 -3.6865 1.681 1.9808 2.2921 -3.3375 -0.9765 2.222 Mean error 0.6782 0.504 0.4425 0.3454 0.415 0.321 0.3071 0.3033 0.4787 0.3932 0.3141 0.4819 0.6001

[0172] Table 5 Statistics of track deviations after 2016 proposed stabilization correction

[0173] flight path 1-2 2-3 3-4 4-5 5-6 6-7 7-8 8-9 9-10 10-11 11-12 12-13 13-14 <![CDATA[RMSE track ]]> 1 1.3167 -0.4535 0.5043 2.0554 -0.9033 3.4754 0.872 2.0473 -3.7385 1.2693 1.0743 -1.5531 -1.0405 2 1.312 -0.2683 0.8635 2.2037 -0.6228 3.2488 1.103 2.1616 -3.8374 1.4959 0.9716 -1.4646 -1.142 3 1.164 -0.4953 0.6771 2.442 -0.8269 3.3776 0.9272 2.3848 -3.5841 1.3776 0.8691 -1.8544 -1.1041 4 1.4015 -0.0425 0.8902 2.2741 -0.7907 3.3712 0.9621 2.1171 -3.5949 1.4902 1.1162 -1.9222 -1.2559 5 1.4997 -0.1786 0.7183 2.1845 -0.9221 3.4156 1.0009 2.3702 -3.7215 1.1174 0.9763 -1.865 -1.1245 6 1.4905 -0.4993 0.7183 2.1042 -0.5905 3.0783 1.0032 2.3464 -3.8685 1.2643 1.1788 -1.79 -1.4026 7 1.0635 -0.4848 0.5246 2.2205 -0.6875 3.2287 0.9931 2.412 -3.6597 1.0257 1.1464 -1.8277 -1.4659 8 1.1161 -0.3958 0.5248 2.4781 -0.6307 3.3091 1.1049 2.414 -3.8832 1.3784 0.8783 -1.4435 -1.4126 9 1.0118 -0.2725 0.5456 2.062 -0.5974 3.4661 0.8834 2.1467 -3.7718 1.301 0.9732 -1.9692 -1.2493 10 1.3037 -0.4364 0.797 2.2354 -0.9664 3.4175 0.894 2.1547 -3.8077 1.4286 0.8733 -1.6444 -1.1021 <![CDATA[ vz ]]> 1.268 -0.3527 0.6764 2.226 -0.7538 3.3388 0.9744 2.2555 -3.7467 1.3148 1.0057 -1.7334 -1.2299 1.885 Mean error 0.5176 0.4673 0.4328 0.4309 0.4363 0.3722 0.2511 0.4255 0.3222 0.4603 0.3495 0.5762 0.4548

[0174] The first phase of point cloud (2009) has the following deviations: Figure 10 As shown in the figure, the flight path deviation of the second phase point cloud (2016) is as follows: Figure 11 As shown in the figure, after the quasi-stabilization adjustment processing of the flight strips in 2009 and 2016, the RMSE track was 2.223 cm and 1.885 cm, respectively, and the RMSE track was reduced by 3.88 cm and 1.93 cm, respectively.

[0175] The method and system of the present invention can improve the relative accuracy of point cloud data between flight strips without flight quality data, making it possible to perform differential calculations of settlement amounts, allowing airborne lidar to play a role in the field of deformation monitoring. In addition, the method has a certain degree of universality and is suitable for point cloud data obtained from most open data sources. There is no need to specifically purchase specific point cloud data, thus reducing research costs.

[0176] The above descriptions are merely a few embodiments of the present application and do not constitute any form of limitation to the present application. Although the present application discloses the preferred embodiments as above, they are not intended to limit the present application. Any technical personnel familiar with the present profession, without departing from the scope of the technical solution of the present application, using the technical content disclosed above to make slight changes or modifications are equivalent to equivalent implementation cases and fall within the scope of the technical solution.

Claims

1. A surface subsidence monitoring method based on laser point cloud data, characterized in that: include: Step 1: Obtain the first phase point cloud data and the second phase point cloud data of the area to be measured in different time periods; Step 2: Determine the flight path deviation of each flight path in the Z coordinate direction in each phase of point cloud data; Step 3: Determine the Z coordinate correction value of each flight strip according to the flight strip deviation, and use the Z coordinate correction value to correct the corresponding flight strip; Step 4: determining the surface settlement of the area to be measured based on the first phase point cloud data and the second phase point cloud data after the corrected flight path; The step 2 specifically includes: Step 21: Obtain the overlapping portion of two adjacent flight strips in each phase of point cloud data, and record each point in the upper flight strip in the overlapping portion as P i , each point in the lower flight path is denoted as P j ; Step 22: Determine each P i and the nearest multiple P j The distance between the fitting planes is recorded as P i The point cloud deviation value of the point; Step 23: Based on multiple P i The point cloud deviation value of the point and the corresponding x-axis coordinate value are used to determine the flight path deviation.

2. The surface subsidence monitoring method based on laser point cloud data according to claim 1, characterized in that: The step 23 specifically includes: According to multiple P i The point cloud deviation value and the corresponding x-axis coordinate value of the point are used to determine the flight path deviation constant using the polynomial fitting method; The flight path deviation of each flight path in the Z coordinate direction in each period of point cloud data is determined according to the flight path deviation constant.

3. The surface subsidence monitoring method based on laser point cloud data according to claim 1, characterized in that: According to the flight path deviation, the Z coordinate correction value of each flight path is determined, including: The Z coordinate correction value of each flight strip is determined according to the first formula, which is: Where V is the deviation error of each flight strip, A is the coefficient matrix, is the Z coordinate correction number of each flight strip, and l is the expected value of the flight strip deviation.

4. The surface subsidence monitoring method based on laser point cloud data according to claim 3 is characterized in that: The Z coordinate correction number Including Z coordinate correction of stable flight zone and Z coordinate correction of unstable flight zone but and Determined according to the second formula, the second formula is: Where M = N 22 -N 21 N 11 - N 12 , is the transpose of A1, - is the sign of the inverse moment, P is the weight matrix, l is the expected value of the flight path deviation, and A1 and A2 are the matrices after the coefficient matrix A is decomposed.

5. The surface subsidence monitoring method based on laser point cloud data according to claim 1, characterized in that: The step 4 specifically includes: The surface settlement of the area to be measured is determined using the M3C2 algorithm based on the first phase point cloud data and the second phase point cloud data after the corrected flight path.

6. The surface subsidence monitoring method based on laser point cloud data according to claim 1, characterized in that: After step 1, the method further includes: Eliminating outliers in the first-phase point cloud data and the second-phase point cloud data; Resampling the first phase point cloud data and the second phase point cloud data after removing outliers using an octree method; Accordingly, step 2 is as follows: Determine the flight path deviation of each flight path in the Z coordinate direction in each phase of the resampled point cloud data.

7. The surface subsidence monitoring method based on laser point cloud data according to claim 6, characterized in that: After resampling the first phase point cloud data and the second phase point cloud data after removing outliers using an octree method, the method further includes: Obtain the point source ID of the midpoint of the first phase point cloud data and the second phase point cloud data after resampling; The first phase point cloud data and the second phase point cloud data are respectively split into flight strips according to the point source ID.

8. The surface subsidence monitoring method based on laser point cloud data according to claim 7, characterized in that: After the first phase point cloud data and the second phase point cloud data are respectively split into flight strips according to the point source ID, the method further includes: The first phase point cloud data and the second phase point cloud data after the flight strip splitting are filtered using a point cloud ground point filtering method to obtain ground points in the first phase point cloud data and ground points in the second phase point cloud data; Accordingly, the step 2 is specifically as follows: Determine the flight path deviation of each flight path in the Z coordinate direction at each ground point.

9. A surface subsidence monitoring system based on laser point cloud data and the method according to any one of claims 1 to 8, characterized in that: include: A data acquisition module, the data acquisition module is used to acquire the first phase point cloud data and the second phase point cloud data of the area to be measured in different time periods; A deviation determination module, the deviation determination module is used to determine the flight path deviation of each flight path in each period of point cloud data in the Z coordinate direction; a flight strip correction module, the flight strip correction module being used to determine a Z coordinate correction number for each flight strip according to the flight strip deviation, and to correct the corresponding flight strip using the Z coordinate correction number; A settlement determination module is used to determine the surface settlement of the area to be measured based on the first phase point cloud data and the second phase point cloud data after the corrected flight path.

Images