Laser point cloud adaptive projection scale-based topographic change detection method

Abstract

The invention discloses a topographic change detection method based on a laser point cloud adaptive projection scale. The method comprises the following steps: 1, obtaining a point cloud of a drainage basin to be detected; 2, obtaining the projection scale of each point in the initial laser point cloud; 3, inputting the projection scale of each point in the initial laser point cloud, and obtaining the M3C2 value of each point in the initial laser point cloud relative to the laser point cloud after the disaster; and 4, performing grid conversion on the M3C2 value of each point in the initial laser point cloud relative to the laser point cloud after the disaster, and obtaining the volume of the accumulation body in the geological disaster area of the drainage basin to be detected based on the multiplication and summation of the grid area and the M3C2 value. The method is simple in step, and the precision of the M3C2 value is improved by introducing the adaptive projection scale, so that the volume of the accumulation body in the geological disaster area of the drainage basin to be detected is obtained, and topographic change detection is realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of watershed terrain change measurement, and particularly relates to a terrain change detection method based on the adaptive projection scale of laser point cloud. Background Technique

[0002] With the continuous development of remote sensing technology, point cloud data, as an important carrier of spatial information, has been widely used in surface disaster monitoring. Point cloud data mainly comes from lidar, which can obtain three-dimensional spatial information with high precision and provides strong data support for terrain change monitoring. However, how to handle the uneven point density, insufficient accuracy in data sparse areas, and point cloud noise and errors in large-scale regions remains a key challenge in the field of terrain change monitoring.

[0003] Currently, existing terrain change detection methods, especially the point cloud to point cloud method. Among them, the M3C2 algorithm can capture the subtle changes in smooth and complex regions of the terrain, enhance the robustness to noise data, and greatly improve the terrain change detection ability. However, when the M3C2 algorithm processes regions with uneven point cloud density, its effectiveness and accuracy decline because the M3C2 algorithm assumes that the point cloud is evenly distributed in space, that is, the projection scale is a fixed value, thus reducing the reliability of the monitoring results.

[0004] Therefore, a terrain change detection method based on the adaptive projection scale of laser point cloud with reasonable design is needed. By introducing the adaptive projection scale, the accuracy of the M3C2 value is improved to obtain the volume of the accumulation body in the geological disaster area of the watershed to be measured, and terrain change detection is realized. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a terrain change detection method based on the adaptive projection scale of laser point cloud for the deficiencies in the above-mentioned prior art. The method steps are simple and reasonably designed. By introducing the adaptive projection scale, the accuracy of the M3C2 value is improved to obtain the volume of the accumulation body in the geological disaster area of the watershed to be measured, and terrain change detection is realized.

[0006] To solve the above technical problem, the technical solution adopted by the present invention is: a terrain change detection method based on the adaptive projection scale of laser point cloud, characterized in that the method includes the following steps:

[0007] Step 1: Acquisition of point cloud of the watershed to be measured:

[0008] Use an unmanned aerial vehicle (UAV) airborne mobile measurement system equipped with a three-dimensional lidar to scan the watershed to be measured to obtain the initial laser point cloud and the post-disaster laser point cloud;

[0009] Step 2: Obtain the projection scale of each point in the initial laser point cloud;

[0010] Step 3: Input the projection scales of the points in the initial lidar point cloud, and obtain the M3C2 values of the points in the initial lidar point cloud relative to the post-disaster lidar point cloud;

[0011] Step 4: Perform raster conversion on the M3C2 values of the points in the initial lidar point cloud relative to the post-disaster lidar point cloud, and sum the products of the raster area and the M3C2 values to obtain the volume of the accumulation body in the geological disaster area of the watershed to be measured.

[0012] The above method for detecting terrain changes based on the adaptive projection scale of lidar point cloud is characterized in that: Step 2, the specific process is as follows:

[0013] Step 201: Use a computer to take any point i in the initial lidar point cloud as the center and r as the search radius to obtain the number of neighboring points N within the search radius of this point; i ; and according to S = πr 2 , obtain the initial projection area S of point i; where π represents the pi;

[0014] Step 202: According to obtain the proportion P of the minimum number of points to the number of neighboring points i ;

[0015] Step 203: Set the projection scale of this point i as d i , then according to obtain the projection area S of point i i ;

[0016] Step 204: According to obtain the projection scale d of point i i ;

[0017] Step 205: According to the method of Step 201 to Step 204, obtain the projection scales of the points in the initial lidar point cloud.

[0018] The above method for detecting terrain changes based on the adaptive projection scale of lidar point cloud is characterized in that: Step 3, the specific process is as follows:

[0019] Step 301: Use a computer to take any point i in the initial lidar point cloud as the center, within the set diameter D range, and use the plane fitting algorithm to obtain the fitting plane of point i and the local normal vector of this fitting plane;

[0020] Step 302: Use a computer to construct a cylinder with a diameter of d i for point i, the center line of the cylinder passes through point i and is along the local normal vector direction; and use a computer to make the diameter d iThe region where the cylinder intersects with the initial laser point cloud is denoted as the intersection point cloud in the q-th period, and the number of points N in the intersection point cloud in the q-th period is obtained. q ; The region where the cylinder with a diameter of d i intersects with the post-disaster laser point cloud is denoted as the intersection point cloud in the (q + 1)-th period, and the number of points N in the intersection point cloud in the (q + 1)-th period is obtained. q+1 ;

[0021] Step 303: Use a computer to judge N q and N q+1 . If N q ≤4 or N q+1 ≤4, then point i is an interference point and is deleted; if 4 < N q and 4 < N q+1 , then execute Step 304;

[0022] Step 304: Use a computer to denote the projection of the line connecting any point in the intersection point cloud in the q-th period and point i on the local normal vector as the projection distance, then obtain the average value of the projection distances of all points in the intersection point cloud in the q-th period, and denote it as the average projection distance L of the intersection point cloud in the q-th period. q ;

[0023] Use a computer to denote the projection of the line connecting any point in the intersection point cloud in the (q + 1)-th period and point i on the local normal vector as the projection distance, then obtain the average value of the projection distances of all points in the intersection point cloud in the (q + 1)-th period, and denote it as the average projection distance L of the intersection point cloud in the (q + 1)-th period. q+1 ;

[0024] Step 305: According to L i = L q - L q+1 , obtain the M3C2 value L i of the initial laser point cloud at point i relative to the post-disaster laser point cloud;

[0025] Step 306: If N q ≥30 and N q+1 ≥30, then according to obtain the absolute value of the uncertainty error LoD; where, σ q,i represents the roughness of the intersection point cloud in the q-th period, σ q+1,i represents the roughness of the intersection point cloud in the (q + 1)-th period, and reg represents the registration error;

[0026] If 4 < N q <30 and 4 < N q+1 <30, then according to obtain the absolute value of the uncertainty error LoD; where, δ q,i represents the standard deviation of the coordinates of each point in the intersection point cloud in the q-th period, δ q+1,iDenote the standard deviation of the coordinates of each point in the intersection point cloud of the (q + 1)-th period;

[0027] Step 307: Compare the M3C2 value L in Step 305 i with the absolute value of the uncertainty error LoD. If L i is greater than LoD, it indicates that the M3C2 value L of the initial laser point cloud at point i relative to the post-disaster laser point cloud i is reliable; otherwise, the M3C2 value L of the initial laser point cloud at point i relative to the post-disaster laser point cloud i is not reliable, and point i is considered an interference point and deleted;

[0028] Step 308: Judge each point of the initial laser point cloud according to the method in Step 301 to Step 307, and obtain the M3C2 value of each point in the initial laser point cloud relative to the post-disaster laser point cloud.

[0029] For the above terrain change detection method based on the adaptive projection scale of the laser point cloud, it is characterized in that: the acquisition of the registration error reg in Step 306 is as follows:

[0030] Step A01: Import the initial laser point cloud and the post-disaster laser point cloud into the CloudCompare software;

[0031] Step A02: Use a computer to select multiple building area point clouds from the initial laser point cloud as the invariant point cloud area of the initial laser point cloud in the CloudCompare software, and select the corresponding building area point cloud from the post-disaster laser point cloud as the invariant point cloud area of the post-disaster laser point cloud;

[0032] Step A03: Use a computer in the CloudCompare software to obtain the M3C2 value of each point in the invariant point cloud area of the initial laser point cloud relative to the invariant point cloud area of the post-disaster laser point cloud, and perform an average value process on it to obtain the registration error reg.

[0033] For the above terrain change detection method based on the adaptive projection scale of the laser point cloud, it is characterized in that: the roughness σ of the q-th period intersection point cloud q,i and the roughness σ of the (q + 1)-th period intersection point cloud q+1,i are obtained as follows:

[0034] Step B01: Obtain the distance a(k) from the k-th point in the q-th period intersection point cloud to the fitting plane in Step 301 and the average value a of the distances from each point in the q-th period intersection point cloud to the fitting plane in Step 301; where k is a positive integer;

[0035] Step B02: According to obtain the roughness σ of the q-th period intersection point cloudq,i ;

[0036] Step B03: Obtain the distance b(e) between the e-th point in the (q + 1)-th intersection point cloud and the fitted plane in Step 301, and the average distance b between each point in the (q + 1)-th intersection point cloud and the fitted plane in Step 301; where e is a positive integer.

[0037] Step B04: According to obtain the roughness σ of the (q + 1)-th intersection point cloud q+1,i .

[0038] The above method for terrain change detection based on the adaptive projection scale of laser point clouds is characterized in that: Step Four, the specific process is as follows:

[0039] Step 401: Use a computer to import the coordinates and M3C2 values of each point in the initial laser point cloud into ArcGIS software, and use the computer to utilize the "Point to Raster" tool in ArcGIS software to set the cell size to 0.1 m, and convert each point data into raster data to obtain a raster map of the watershed to be measured; where the raster map of the watershed to be measured is in the.tif format.

[0040] Step 402: Conduct on-site investigation of the geological disaster area, and use a GNSS-RTK device to obtain the coordinates of each point on the boundary line of the accumulation body in the geological disaster area.

[0041] Step 403: Input the coordinates of each point on the boundary line of the accumulation body in the geological disaster area into ArcGIS software to obtain the area of the accumulation body surface.

[0042] According to the area of the accumulation body surface, use a computer to utilize the "Extract by Mask" tool in ArcGIS software to obtain the raster area of the accumulation body corresponding to the area of the accumulation body surface from the raster map of the watershed to be measured.

[0043] Step 404: Multiply the area of each grid cell in the raster area of the accumulation body by the absolute value of the M3C2 value and sum them up to obtain the volume of the accumulation body in the geological disaster area.

[0044] The present invention has the following advantages compared with the prior art:

[0045] 1. The method of the present invention has simple steps, reasonable design and convenient implementation, and high accuracy.

[0046] 2. The present invention uses an unmanned aerial vehicle (UAV) airborne mobile measurement system equipped with a three-dimensional lidar to scan the watershed area to be measured, and obtains the initial laser point cloud and the post-disaster laser point cloud, which is convenient for subsequent terrain change detection based on the initial laser point cloud and the post-disaster laser point cloud.

[0047] 3. The present invention obtains the projection scales of each point in the initial laser point cloud. By introducing an adaptive projection scale, the accuracy of subsequent M3C2 values is improved, thus adapting to the sparse or dense phenomenon of the point cloud distribution.

[0048] 4. The present invention performs raster conversion on the M3C2 values of each point in the initial laser point cloud relative to the post-disaster laser point cloud, and sums the products of the raster area and the M3C2 values to obtain the volume of the accumulation body in the geological disaster area of the basin to be measured, realizing terrain change detection.

[0049] In summary, the method steps of the present invention are simple and reasonably designed. By introducing an adaptive projection scale, the accuracy of M3C2 values is improved to obtain the volume of the accumulation body in the geological disaster area of the basin to be measured and realize terrain change detection.

[0050] The technical solution of the present invention will be further described in detail below with reference to the drawings and embodiments. Description of the Drawings

[0051] Figure 1 is the flowchart of the method of the present invention.

[0052] Figure 2 is the schematic diagram of each point of the boundary line of the accumulation body in the geological disaster area surveyed on the spot of the present invention. Detailed Embodiments

[0053] As Figure 1 shown, a terrain change detection method based on the adaptive projection scale of laser point cloud includes the following steps:

[0054] Step 1: Acquisition of point cloud of the basin to be measured:

[0055] Use an unmanned aerial vehicle (UAV) airborne mobile measurement system equipped with a three-dimensional lidar to scan the basin to be measured to obtain the initial laser point cloud and the post-disaster laser point cloud;

[0056] Step 2: Obtain the projection scale of each point in the initial laser point cloud;

[0057] Step 3: Input the projection scales of each point in the initial laser point cloud to obtain the M3C2 values of each point in the initial laser point cloud relative to the post-disaster laser point cloud;

[0058] Step 4: Perform raster conversion on the M3C2 values of each point in the initial laser point cloud relative to the post-disaster laser point cloud, and sum the products of the raster area and the M3C2 values to obtain the volume of the accumulation body in the geological disaster area of the basin to be measured.

[0059] In this embodiment, the specific process of Step 2 is as follows:

[0060] Step 201: Using a computer, in the initial laser point cloud, with any point i as the center and r as the search radius, obtain the number of neighboring points N within the search radius of this point; and according to S = πr², obtain the initial projected area S of point i; where π represents the pi. i ; and according to S = πr 2 , obtain the initial projected area S of point i; where π represents the pi.

[0061] Step 202: According to obtain the ratio P of the minimum number of points to the number of neighboring points i ;

[0062] Step 203: Set the projected scale of this point i as d i , then according to obtain the projected area S of point i i ;

[0063] Step 204: According to obtain the projected scale d of point i i ;

[0064] Step 205: According to the methods of Step 201 to Step 204, obtain the projected scales of each point in the initial laser point cloud.

[0065] In this embodiment, Step 3 is specifically as follows:

[0066] Step 301: Using a computer, in the initial laser point cloud, with any point i as the center, within the set diameter D range, adopt the plane fitting algorithm to obtain the fitting plane of point i and the local normal vector of this fitting plane;

[0067] Step 302: Using a computer, construct a cylinder with a diameter of d i for point i, the center line of the cylinder passes through point i and is along the local normal vector direction; and use the computer to record the area where the cylinder with a diameter of d i intersects with the initial laser point cloud as the q-th phase intersection point cloud, and obtain the number of points N of the q-th phase intersection point cloud q ; record the area where the cylinder with a diameter of d i intersects with the post-disaster laser point cloud as the (q + 1)-th phase intersection point cloud, and obtain the number of points N of the (q + 1)-th phase intersection point cloud q+1 ;

[0068] Step 303: Using a computer, judge N q and N q+1 , if N q ≤4 or N q+1 ≤4, then point i is an interference point and is deleted; if 4 < N q and 4 < N q+1 , then execute Step 304;

[0069] Step 304: Use a computer to record the projection of the line connecting any point in the q-th intersection point cloud to point i on the local normal vector as the projection distance, and then obtain the average value of the projection distances of all points in the q-th intersection point cloud, which is denoted as the projection distance mean L of the q-th intersection point cloud. q ;

[0070] Use a computer to record the projection of the line connecting any point in the (q + 1)-th intersection point cloud to point i on the local normal vector as the projection distance, and then obtain the average value of the projection distances of all points in the (q + 1)-th intersection point cloud, which is denoted as the projection distance mean L of the (q + 1)-th intersection point cloud. q+1 ;

[0071] Step 305: According to L i = L q - L q+1 , obtain the M3C2 value L of the initial laser point cloud at point i relative to the post-disaster laser point cloud. i ;

[0072] Step 306: If N q ≥ 30 and N q+1 ≥ 30, then according to , obtain the absolute value of the uncertainty error LoD; where, σ q,i represents the roughness of the q-th intersection point cloud, σ q+1,i represents the roughness of the (q + 1)-th intersection point cloud, and reg represents the registration error.

[0073] If 4 < N q < 30 and 4 < N q+1 < 30, then according to , obtain the absolute value of the uncertainty error LoD; where, δ q,i represents the standard deviation of the coordinates of each point in the q-th intersection point cloud, and δ q+1,i represents the standard deviation of the coordinates of each point in the (q + 1)-th intersection point cloud.

[0074] Step 307: Compare the M3C2 value L i in Step 305 with the absolute value of the uncertainty error LoD. If L i is greater than LoD, it means that the M3C2 value L i of the initial laser point cloud at point i relative to the post-disaster laser point cloud is credible; otherwise, the M3C2 value L i of the initial laser point cloud at point i relative to the post-disaster laser point cloud is not credible, and point i is considered an interference point and deleted.

[0075] Step 308: According to the method in Step 301 to Step 307, judge each point of the initial laser point cloud to obtain the M3C2 value of each point in the initial laser point cloud relative to the post-disaster laser point cloud.

[0076] In this embodiment, the process of obtaining the registration error reg in step 306 is as follows:

[0077] Step A01: Import the initial laser point cloud and the post-disaster laser point cloud into the CloudCompare software;

[0078] Step A02: Use a computer to select multiple building area point clouds from the initial laser point cloud in the CloudCompare software as the invariant point cloud area of the initial laser point cloud, and select the corresponding building area point cloud from the post-disaster laser point cloud as the invariant point cloud area of the post-disaster laser point cloud;

[0079] Step A03: Use a computer in the CloudCompare software to obtain the M3C2 values of each point in the invariant point cloud area of the initial laser point cloud relative to the invariant point cloud area of the post-disaster laser point cloud, and perform an average processing on them to obtain the registration error reg.

[0080] In this embodiment, the roughness σ of the q-th intersection point cloud q,i and the roughness σ of the (q + 1)-th intersection point cloud q+1,i are obtained as follows:

[0081] Step B01: Obtain the distance a(k) from the k-th point in the q-th intersection point cloud to the fitting plane in step 301 and the average distance a of each point in the q-th intersection point cloud to the fitting plane in step 301; where k is a positive integer;

[0082] Step B02: According to obtain the roughness σ of the q-th intersection point cloud q,i ;

[0083] Step B03: Obtain the distance b(e) from the e-th point in the (q + 1)-th intersection point cloud to the fitting plane in step 301 and the average distance b of each point in the (q + 1)-th intersection point cloud to the fitting plane in step 301; where e is a positive integer;

[0084] Step B04: According to obtain the roughness σ of the (q + 1)-th intersection point cloud q+1,i ;

[0085] As Figure 2 shown, in this embodiment, step four is as follows:

[0086] Step 401: Use a computer to import the coordinates of each point and the M3C2 values in the initial laser point cloud into ArcGIS software, and use the computer to utilize the "Point to Raster" tool in ArcGIS software to set the cell size to 0.1 m, and convert the data of each point into raster data to obtain a raster map of the basin to be measured; among them, the raster map of the basin to be measured is in the.tif format;

[0087] Step 402: Conduct on-site inspections of the geological disaster area, and use a GNSS-RTK device to obtain the coordinates of each point on the boundary line of the accumulation body in the geological disaster area;

[0088] Step 403: Input the coordinates of each point on the boundary line of the accumulation body in the geological disaster area into ArcGIS software to obtain the area of the accumulation body surface;

[0089] According to the area of the accumulation body surface, use a computer to utilize the "Extract by Mask" tool in ArcGIS software to obtain the raster area of the accumulation body corresponding to the area of the accumulation body surface from the raster map of the basin to be measured;

[0090] Step 404: Multiply the area of each grid cell in the raster area of the accumulation body by the absolute value of the M3C2 value and sum them to obtain the volume of the accumulation body in the geological disaster area.

[0091] In this embodiment, the accumulation body refers to the area where sedimentation occurs in the geological disaster area.

[0092] In this embodiment, the acquisition of the standard deviation δ of the coordinates of each point in the intersecting point cloud is specifically as follows:

[0093] where x n , y n , z n represent the X coordinate, Y coordinate, and Z coordinate of the nth point in the intersecting point cloud, represent the average values of the X coordinates, Y coordinates, and Z coordinates of N points in the intersecting point cloud, 1 ≤ n ≤ N, and N is the total number of points in the intersecting point cloud.

[0094] In this embodiment, use the SZT-R250 UAV airborne mobile measurement system equipped with a RIEGL VUX-1UAV 3D lidar to collect scanning data along the set flight path; among them, the flight altitude is 70 m higher than the terrain altitude, the flight speed is 7 km / h, the scanning strip width is 60 m, the scanning angle is 360°, and the scanning pulse rate is 100 kHz.

[0095] In this embodiment, the set flight path is planned by using the "DOUBLE GRID For 3D Models" function provided by pix4D software, and the overall flight path is in the shape of a "well" character.

[0096] In this embodiment, on-site field surveys need to be carried out before setting the flight route.

[0097] In this embodiment, in step 1, the UAV-borne mobile measurement system is equipped with a 3D lidar to scan the basin to be measured. The specific process is as follows:

[0098] Step 101: Obtain the first flight of lidar point clouds,..., the f-th flight of lidar point clouds,..., the F-th flight of lidar point clouds; where f and F are both positive integers, 1 ≤ f ≤ F, and F ≥ 5;

[0099] Step 102: Use a computer to utilize the "Merge" tool in CloudCompare software to merge the lidar point clouds from the first flight to the F-th flight to obtain the merged initial lidar point cloud;

[0100] Step 103: Use a computer to utilize the "Segment" tool in CloudCompare software to crop the merged initial lidar point cloud according to the results of the on-site field surveys of the basin to be measured to obtain the point cloud of the basin to be measured;

[0101] Step 104: Use a computer to perform a first filtering on the point cloud of the basin to be measured using the MCC point cloud filtering method to obtain the initial lidar point cloud after the first filtering;

[0102] Step 105: Use a computer to perform a second filtering on the initial lidar point cloud after the first filtering using TerraSolid software to obtain the initial lidar point cloud after the second filtering, and denote it as the initial lidar point cloud;

[0103] Step 106: After a geological disaster occurs in the basin to be measured, according to the methods described in steps 101 to 105, obtain the post-disaster lidar point cloud.

[0104] In this embodiment, the coordinates of each point in the initial lidar point cloud and the post-disaster lidar point cloud are in the CGCS2000 coordinate system. The coordinates of the boundary points obtained by the GNSS-RTK device are also in the CGCS2000 coordinate system.

[0105] In this embodiment, the search radius r takes a value of 1 m or 0.5 m; where when r is 1 m, it is suitable for flat terrain areas; when r is 0.5 m, it is suitable for hilly terrain areas.

[0106] In this embodiment, the set diameter D in step 301 takes a value of 2d i 。

[0107] In summary, the design of the present invention is reasonable. By introducing an adaptive projection scale, the accuracy of the M3C2 value is improved, so as to obtain the volume of the accumulation body in the geological disaster area of the basin to be measured and realize terrain change detection.

[0108] The above are only the preferred embodiments of the present invention and do not impose any limitations on the present invention. Any simple modifications, changes, and equivalent structural changes made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A terrain change detection method based on adaptive projection scale of laser point cloud, characterized in that: The method comprises the following steps: Step 1: Obtain the point cloud of the watershed to be measured: The UAV-mounted mobile measurement system equipped with a three-dimensional laser radar is used to scan the watershed to obtain the initial laser point cloud and the post-disaster laser point cloud; Step 2: Obtain the projection scale of each point in the initial laser point cloud; Step 3: Input the projection scale of each point in the initial laser point cloud to obtain the M3C2 value of each point in the initial laser point cloud relative to the laser point cloud after the disaster; Step 4: Perform raster conversion on each point in the initial laser point cloud relative to the M3C2 value of the laser point cloud after the disaster, and obtain the volume of the accumulation body in the geological disaster area of ​​the measured basin based on the multiplication and summation of the grid area and the M3C2 value.

2. A terrain change detection method based on laser point cloud adaptive projection scale according to claim 1, characterized in that: Step 2: The specific process is as follows: Step 201: Use a computer to obtain the number of neighboring points N within the search radius of any point i in the initial laser point cloud, with r as the search radius. i ; And according to S = πr 2 , get the initial projection area S of point i; where π represents pi; Step 202: According to Get the ratio P of the minimum number of points to the number of adjacent points i ; Step 203: Set the projection scale of the point i to d i , then according to Get the projection area S of point i i ; Step 204: According to Get the projection scale d of point i i ; Step 205: According to the method of steps 201 to 204, the projection scale of each point in the initial laser point cloud is obtained.

3. A terrain change detection method based on laser point cloud adaptive projection scale according to claim 1, characterized in that: Step 3: The specific process is as follows: Step 301: Using a computer, in the initial laser point cloud, taking any point i as the center and within a set diameter D, using a plane fitting algorithm, obtain a fitting plane of point i and a local normal vector of the fitting plane; Step 302: Use a computer to construct a diameter d for point i. i The center line of the cylinder passes through point i and is along the direction of the local normal vector; and a computer is used to calculate the diameter d i The area where the cylinder intersects with the initial laser point cloud is recorded as the qth period intersection point cloud, and the number of points N of the qth period intersection point cloud is obtained q ; Set the diameter to d i The area where the cylinder intersects with the laser point cloud after the disaster is recorded as the q+1th period intersection point cloud, and the number of points N of the q+1th period intersection point cloud is obtained q+1 ; Step 303: Use a computer to q and N q+1 Make a judgment, if N q ≤4 or N q+1 ≤4, then point i is an interference point and is deleted; if 4 <N q And 4 <N q+1 , then execute step 304; Step 304: Use a computer to record the projection of the line connecting any point in the qth period intersection point cloud and point i on the local normal vector as the projection distance, and then obtain the average value of the projection distances of all points in the qth period intersection point cloud, and record it as the projection distance mean value L of the qth period intersection point cloud. q ; The projection of the line connecting any point in the q+1th intersecting point cloud and point i on the local normal vector is recorded as the projection distance by a computer, and then the average value of the projection distances of all points in the q+1th intersecting point cloud is obtained and recorded as the projection distance mean value L of the q+1th intersecting point cloud. q+1 ; Step 305: According to L i =L q -L q+1 , and obtain the M3C2 value L of the laser point cloud at the initial laser point cloud point i relative to the laser point cloud after the disaster i ; Step 306: If N q ≥30 and N q+1 ≥30, then according to Get the absolute value of uncertainty error LoD; where σ q,i represents the roughness of the qth intersection point cloud, σ q+1,i represents the roughness of the intersection point cloud at the q+1th period, and reg represents the registration error; If 4 <N q <30 and 4 <N q+1 <30, then according to Get the absolute value of uncertainty error LoD; where δ q,i represents the standard deviation of the coordinates of each point in the qth intersecting point cloud, δ q+1,i Represents the standard deviation of the coordinates of each point in the intersection point cloud of the q+1th period; Step 307: Set the M3C2 value in step 305 to L i Compared with the absolute value of uncertainty error LoD, if L i If it is greater than LoD, it means that the M3C2 value L of the initial laser point cloud point i relative to the laser point cloud after the disaster i Credible; otherwise, the M3C2 value L of the initial laser point cloud point i relative to the laser point cloud after the disaster i If it is untrustworthy, point i is considered as an interference point and deleted; Step 308: According to the method of steps 301 to 307, each point of the initial laser point cloud is judged to obtain the M3C2 value of each point in the initial laser point cloud relative to the laser point cloud after the disaster.

4. A terrain change detection method based on laser point cloud adaptive projection scale according to claim 3, characterized in that: In step 306, the registration error reg is obtained, and the specific process is as follows: Step A01, importing the initial laser point cloud and the post-disaster laser point cloud into CloudCompare software; Step A02, using a computer to use CloudCompare software to select multiple building area point clouds from the initial laser point cloud as the initial laser point cloud unchanged point cloud area, and select corresponding building area point clouds from the post-disaster laser point cloud as the post-disaster laser point cloud unchanged point cloud area; Step A03, using a computer to use CloudCompare software to obtain the M3C2 value of each point in the unchanged point cloud area of ​​the initial laser point cloud relative to the unchanged point cloud area of ​​the laser point cloud after the disaster, and average the values ​​to obtain the registration error reg.

5. A terrain change detection method based on laser point cloud adaptive projection scale according to claim 3, characterized in that: The roughness σ of the qth intersection point cloud q,i The roughness σ of the point cloud intersecting with the q+1th period q+1,i The specific process of obtaining is as follows: Step B01, obtaining the distance a(k) between the kth point in the qth period intersection point cloud and the fitting plane in step 301 and the average distance a between each point in the qth period intersection point cloud and the fitting plane in step 301; wherein k is a positive integer; Step B02: Get the roughness σ of the qth intersection point cloud q,i ; Step B03, obtaining the distance b(e) between the e-th point in the q+1-th intersecting point cloud and the fitting plane in step 301 and the average value b of the distances between each point in the q+1-th intersecting point cloud and the fitting plane in step 301; wherein e is a positive integer; Step B04: Get the roughness σ of the q+1th intersection point cloud q+1,i .

6. A terrain change detection method based on laser point cloud adaptive projection scale according to claim 1, characterized in that: Step 4: The specific process is as follows: Step 401: using a computer to import the coordinates and M3C2 values ​​of each point in the initial laser point cloud into ArcGIS software, and using the "point to raster" tool in ArcGIS software, using a computer to set the pixel size to 0.1m, converting each point data into raster data to obtain a raster map of the watershed to be measured; wherein the raster map of the watershed to be measured is in .tif format; Step 402: Conduct an on-site survey of the geological disaster area, and use GNSS-RTK equipment to obtain the coordinates of each point on the boundary line of the accumulation body in the geological disaster area; Step 403, input the coordinates of each point of the boundary line of the accumulation body in the geological disaster area into ArcGIS software to obtain the accumulation body surface area; According to the surface area of ​​the accumulation body, the "Extract by Mask" tool in ArcGIS software is used to obtain the accumulation body grid area corresponding to the surface area of ​​the accumulation body from the raster map of the watershed to be tested; Step 404: multiply the area of ​​each grid unit in the accumulation body grid area by the absolute value of the M3C2 value and sum them up to obtain the volume of the accumulation body in the geological disaster area.

Images