Ground object three-dimensional change detection method, electronic device and storage medium

Abstract

The present invention discloses a method for detecting three-dimensional changes in ground objects, an electronic device and a storage medium, wherein the detection method comprises accurately registering two time-phase point clouds, and setting a reference time-phase point cloud and a time-phase point cloud to be detected; calculating the normal vectors of the time-phase point cloud to be detected and the reference time-phase point cloud; calculating the local Euclidean distance and the local density; improving the local Euclidean distance according to the local density to obtain the density-adaptive local Euclidean distance; extracting the points to be detected whose density-adaptive local Euclidean distance in the time-phase point cloud to be detected is greater than the distance threshold to the candidate change area; eliminating the misdetected points in the candidate change area caused by modeling deformation; clustering the points in the candidate change area, calculating the area corresponding to the clustering cluster, and further eliminating the misdetected points according to the area and the area threshold. The present invention has better noise resistance and better detection accuracy.

Description

Technical Field

[0001] The present invention belongs to the field of surveying and mapping science and technology, and in particular relates to a method for detecting three-dimensional changes of ground objects based on point cloud density adaptive local Euclidean distance, an electronic device and a storage medium. Background Art

[0002] Because point clouds contain elevation information, they can reflect the geometric properties of objects. Geometric changes in objects are often less likely to be misjudged and more pronounced (e.g., new construction or demolition). Therefore, point cloud-based 3D change detection is a good way to capture this information. Furthermore, with the recent development of sensor technology and drones, the accuracy and quality of point clouds are increasing, the information they contain is becoming richer, and the barriers to acquisition and costs are significantly decreasing. Consequently, an increasing number of researchers are conducting research on 3D change detection.

[0003] According to the detection unit, three-dimensional change detection can be mainly divided into two types of detection methods: point-by-point and object-based. Point-by-point 3D change detection usually uses point-by-point height difference, Euclidean distance or graph cut method to extract candidate change areas. On this basis, geometric structure, texture or color information is used for post-refinement processing, such as Murakami H, Nakagawa K, Hasegawa H, et al. Change detection of buildings using an airborne laserscanner [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 1999, 54 (2-3): 148-152., Jung F. Detecting building changes from multitemporal aerialstereopairs [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2004, 58 (3-4): 187-201., Peng Daifeng, Zhang Yongjun, Xiong Xiaodong. 3D change detection of buildings using LiDAR point cloud and aerial imagery [J]. Journal of Wuhan University (Information Science Edition), 2015, 40 (4): 462-468., and Du S, Zhang Y, Qin R, et al. Building change detection using old aerial images and new LiDAR data[J].Remote Sensing,2016,8(12):1030..

[0004] Object-based 3D change detection methods usually require classifying data at different time phases and then comparing the classification results of point clouds to extract the changed areas, such as Matikainen L, J,Ahokas E,etal.Automatic detection of buildings and changes in buildings for updating ofmaps[J].Remote Sensing,2010,2(5):1217-1248.,Pang S,Hu X,Wang Z,et al.Object-based analysis of airborne LiDAR data for building change detection[J].RemoteSensing,2014,6(11):10733-10749., and Qin R, Huang

[0005] Point-by-point 3D change detection methods have high requirements for point cloud quality and registration. Furthermore, seasonal changes in vegetation can lead to false detections, limiting these methods. Compared to point-by-point 3D change detection, object-based 3D change detection has relatively low requirements for point cloud quality and registration. These methods separate the classification and change detection processes, causing the final change detection results to be affected by previous classification errors, leading to error accumulation and a dependence of detection accuracy on the quality of point cloud classification.

[0006] In general, most current three-dimensional change detection methods can achieve high change detection accuracy when the point cloud geometric accuracy is high, the registration error is small, and the objects in the scene change significantly. However, if there is object deformation in the point cloud or a certain registration error, it will lead to a large number of false detections and missed detections of objects, and the change detection accuracy is difficult to guarantee. Summary of the Invention

[0007] The purpose of the present invention is to provide a method, electronic device and storage medium for detecting three-dimensional changes in ground objects, so as to solve the problems of false detection and missed detection when the geometric accuracy of point clouds is not high and the registration error is large.

[0008] The present invention solves the above technical problems through the following technical solutions: A method for detecting three-dimensional changes in ground objects, comprising the following steps:

[0009] Step 1: Read the point clouds of two time phases, perform precise alignment on them, and use one of the two time phase point clouds after precise alignment as the reference time phase point cloud, and the other as the time phase point cloud to be detected;

[0010] Step 2: Calculate the normal vectors of the time phase point cloud to be detected and the reference time phase point cloud;

[0011] Step 3: Calculate each point to be detected in the time phase point cloud to be detected The angle θ between the normal vector and the z direction;

[0012] If the angle θ is less than the angle threshold T angle , then search for the point to be detected Find k neighboring points on the xy plane in the reference time phase point cloud, and cluster the k neighboring points to obtain the first clustering result. According to the first clustering result, select the appropriate cluster to calculate the local Euclidean distance and local density; otherwise, search for the point to be detected Find k neighboring points in the three-dimensional space of the reference phase point cloud and calculate the local Euclidean distance and local density;

[0013] Step 4: improving the local Euclidean distance according to the local density to obtain a density-adaptive local Euclidean distance;

[0014] Step 5: The density adaptive local Euclidean distance in the phase point cloud to be detected is greater than the distance threshold T d The points to be detected are extracted into candidate change areas;

[0015] Step 6: Search for candidate change points q in the candidate change region i k adjacent points on the xy plane in the phase point cloud to be detected are clustered to obtain a second clustering result, and the candidate change point q is clustered according to the second clustering result. i Performing a process of moving out of or not moving out of the candidate change area;

[0016] Step 7: Cluster the points in the candidate change area processed in step 6, and calculate the area S corresponding to each cluster according to the average density of each cluster; if the area S is greater than the area threshold T area , then the cluster is considered to have changed; otherwise, the cluster is considered to be a false detection.

[0017] Preferably, in step 1, the ICP algorithm is used to perform precise registration on the point clouds of the two time phases.

[0018] Furthermore, in step 2, the specific implementation process of calculating the normal vectors of the time phase point cloud to be detected and the reference time phase point cloud is:

[0019] Establishing a three-dimensional KD tree of the time phase point cloud to be detected and a three-dimensional KD tree of the reference time phase point cloud;

[0020] The normal vectors of the to-be-detected time phase point cloud and the reference time phase point cloud are respectively calculated according to the three-dimensional KD tree.

[0021] Furthermore, in step 3, the angle θ is calculated as follows:

[0022]

[0023] Among them, n x Point to be tested The value of the normal vector in the x direction, n y Point to be tested The value of the normal vector in the y direction.

[0024] Furthermore, in step 3, the specific implementation process of selecting a suitable cluster to calculate the local Euclidean distance and local density is as follows:

[0025] Condition A: If the number of clusters in the first clustering result is 1, then the k neighboring points are placed into the parameter calculation point cluster;

[0026] Condition B: If the number of clusters in the first clustering result is 2, and the angle between the normal vector of the point with the largest elevation value among the k neighboring points and the z direction is less than the angle threshold T angle , then the farthest away from the detection point The clusters are put into the parameter calculation point clusters;

[0027] If conditions A and B are not met, the nearest neighbor detection point The clusters are put into the parameter calculation point clusters;

[0028] Calculate the local Euclidean distance of the points in the point cluster according to the parameters;

[0029] Calculate the parameters to calculate the center of the cluster and search for k in the reference phase point cloud d The local density is calculated based on the neighboring points.

[0030] Preferably, the calculation formula of the local Euclidean distance is:

[0031]

[0032] Among them, d Loc is the local Euclidean distance, k r Calculate the number of neighboring points in the cluster for the parameter, Point to be tested At the jth neighbor point in the reference phase point cloud, (x i ,y i ,z i ) is the point to be detected The three-dimensional coordinates, (x j ,y j ,z j ) is a point The three-dimensional coordinates of .

[0033] Preferably, the calculation formula of the local density is:

[0034]

[0035] Among them, ρ Loc is the local density, k d is the number of neighboring points in the reference phase point cloud, r max Calculate the center of the cluster to k for the parameter d The maximum distance between neighboring points.

[0036] Furthermore, in step 4, the calculation formula of the density-adaptive local Euclidean distance is:

[0037]

[0038] in, is the density-adaptive local Euclidean distance, d Loc is the local Euclidean distance, ρ Loc is the local density.

[0039] Furthermore, in step 6, the candidate change point q is clustered according to the second clustering result. i The specific implementation process of moving out or not moving out of the candidate change area is as follows:

[0040] If the number of clusters in the second clustering result is 1, then the candidate change point q i It does not belong to the area where the ground object modeling is poor, and the candidate change point q i Carry out removal processing;

[0041] If the number of clusters in the second clustering result is 2, 3 or 4, then the candidate change point q i Belong to the area where the ground object modeling is not good, if the candidate change point q i The first top point of exists and has not changed, then the candidate change point q i Move out of the candidate change area; otherwise, do not select the candidate change point q i Processing; Among them, the candidate change point q iThe first top point is the point where the angle θ between the normal vector of the k neighboring points and the z direction is less than the angle threshold T angle And the point with the largest elevation value;

[0042] If the number of clusters in the second clustering result is greater than 4, the candidate change point q i Belong to the area where the ground object modeling is not good, if the candidate change point q i If the second top point of the candidate change point q changes, i Process; otherwise, the candidate change point q i Move out of the candidate change area; among them, the candidate change point q i The second top point of is the point with the largest elevation value among the k neighboring points.

[0043] Furthermore, in step 7, the average density of each cluster is calculated as follows:

[0044]

[0045] in, is the average density of the current cluster, k c is the number of points in the current cluster, (x j ,y j ,z j ) is the three-dimensional coordinate of the j-th point in the current cluster, (x n ,y n ,z n ) is the three-dimensional coordinate of the nearest point corresponding to the j-th point in the current cluster.

[0046] Furthermore, in step 7, the calculation formula of area S is:

[0047]

[0048] Among them, S is the area of ​​the current cluster, is the average density of the current cluster, k c is the number of points in the current cluster.

[0049] The present invention also provides an electronic device comprising a memory and a processor, wherein the memory stores a computer program that can be run on the processor, and when the processor runs the computer program, the steps of the above-mentioned method for detecting three-dimensional changes in ground objects are executed.

[0050] The present invention also provides a computer-readable storage medium, which is a non-volatile storage medium or a non-transient storage medium, on which a computer program is stored. When the computer program is run by a processor, the steps of the three-dimensional change detection method of the above-mentioned ground object are executed.

[0051] Beneficial effects

[0052] Compared with the prior art, the advantages of the present invention are:

[0053] The present invention provides a method, electronic device, and storage medium for detecting three-dimensional changes in land objects. Based on the Euclidean distance between a point to be detected and a single neighboring point in three-dimensional space, a local Euclidean distance is introduced between the point to be detected and k neighboring points on the xy plane and in three-dimensional space. This not only reduces false detection of disappearing land objects at the edge, but also greatly improves noise resistance and adaptability. Local density is used to improve the local Euclidean distance, reducing the influence of inconsistent point cloud density and registration errors that cause the calculated distance to be greater than the actual change distance. Neighboring points on the xy plane of a candidate change area are clustered, and whether the current candidate change point belongs to a poorly modeled area of ​​the land object is determined based on the clustering results, thereby reducing false detection caused by point cloud deformation. Local density is used to calculate the area of ​​clusters of candidate change areas, and candidate change areas with smaller areas are eliminated, thereby improving the quality of change detection.

[0054] The method described in the present invention has good anti-noise performance and better change detection effect, and can greatly reduce the false detection caused by deformation of ground object modeling; it can extract the changed area in the point cloud more completely and with higher precision when the geometric accuracy of the point cloud is not high and there is a certain alignment error. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only one embodiment of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0056] Figure 1 is a flow chart of a method for detecting three-dimensional changes in ground objects according to an embodiment of the present invention;

[0057] Figure 2 is the detection data used in the embodiments of the present invention, wherein Figure (a) is the building point cloud in June 2020, and Figure (b) is the building point cloud in September 2020;

[0058] Figure 3 This is the three-dimensional change detection result in the embodiment of the present invention, where Figure (a) is based on Figure 2 (b) is the reference time point cloud, Figure 2 (a) is the change detection result obtained by the point cloud of the time phase to be detected; Figure (b) is based on Figure 2 (a) is the reference time phase point cloud, and (b) is the change detection result obtained from the time phase point cloud to be detected;

[0059] Figure 4 It is the manually labeled accuracy verification data used in the embodiments of the present invention. DETAILED DESCRIPTION

[0060] The following is a clear and complete description of the technical solutions of the present invention in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts are within the scope of protection of the present invention.

[0061] The following specific embodiments are used to describe the technical solution of the present application in detail. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0062] like Figure 1 As shown, a method for detecting three-dimensional changes in ground objects provided by an embodiment of the present invention includes the following steps:

[0063] Step 1: Point cloud reading and precise registration

[0064] Read the point clouds of the two time phases, perform precise registration on the point clouds of the two time phases, and use one of the point clouds of the two time phases after precise registration as the reference time phase point cloud, and the other as the time phase point cloud to be detected. Figure 2 As shown, Figure (a) shows the building point cloud in June 2020, and Figure (b) shows the building point cloud in September 2020. Figure (a) can be used as the reference time phase point cloud, and Figure (b) can be used as the time phase point cloud to be detected, or Figure (b) can be used as the reference time phase point cloud, and Figure (a) can be used as the time phase point cloud to be detected.

[0065] In this embodiment, the ICP algorithm is used to precisely register point clouds from two time phases. The ICP algorithm is simple and easy to implement, and with good initial values, it converges quickly and accurately. Using the ICP algorithm to precisely register point clouds from two time phases is known in the art. For example, see Besl PJ, McKay ND. Method for registration of 3-D shapes [C] / / Sensor fusion IV: control paradigms and data structures. Spie, 1992, 1611: 586-606.

[0066] Step 2: Calculate the normal vectors of the detection phase point cloud and the reference phase point cloud. The specific implementation process is as follows:

[0067] Step 2.1: Establish a 3D KD tree of the time phase point cloud to be detected and a 3D KD tree of the reference time phase point cloud;

[0068] Step 2.2: Calculate the normal vector of the point cloud of the phase to be detected based on the 3D KD tree of the point cloud of the phase to be detected, and calculate the normal vector of the point cloud of the reference phase based on the 3D KD tree of the point cloud of the reference phase.

[0069] Calculating the normal vector of a temporal point cloud using a 3D KD tree is an existing technique. For example, see Holz D, Holzer S, Rusu R B, et al. Real-time plane segmentation using RGB-D cameras [C] / / Robot Soccer World Cup. Springer, Berlin, Heidelberg, 2011: 306-317.

[0070] Step 3: Calculation of local Euclidean distance and local density

[0071] Calculate each point to be detected in the point cloud of the time phase to be detected The angle θ between the normal vector and the z direction; if the angle θ is less than the angle threshold T angle , then search for the point to be detected In the reference time phase point cloud, k adjacent points on the xy plane are clustered and the k adjacent points are clustered to obtain the first clustering result. According to the first clustering result, the appropriate cluster is selected to calculate the local Euclidean distance and local density; otherwise, the point to be detected is searched. Find k neighboring points in the three-dimensional space of the reference phase point cloud and calculate the local Euclidean distance and local density.

[0072] In this embodiment, the x and y coordinates of the time phase point cloud to be detected and the reference time phase point cloud are obtained based on WGS84 using UTM projection, and the z coordinate is the elevation of the WGS84 coordinate system.

[0073] In this embodiment, the point to be detected The angle θ between the normal vector and the z direction is calculated as follows:

[0074]

[0075] Among them, n x Point to be tested The value of the normal vector in the x direction, n y Point to be tested The value of the normal vector in the y direction. Angle threshold T angle It can be set according to the needs. When the point cloud geometric accuracy is high and the registration is good, if you are interested in targets with smaller changes, the angle threshold T angle It can be set smaller, otherwise it can be set larger. In this embodiment, the angle threshold T angleSet to 20°.

[0076] In order to improve the search speed, a two-dimensional KD tree of the phase point cloud to be detected and the reference phase point cloud is established. When searching among k neighboring points, the two-dimensional KD tree based on the reference phase point cloud can quickly search for the point to be detected. k neighboring points on the xy plane of the reference point cloud. Since neighboring points found in the xy plane of the reference point cloud have multiple elevation levels, Euclidean clustering is required. However, neighboring points found in the three-dimensional space of the reference point cloud have only one elevation level and do not require Euclidean clustering. In this embodiment, k is 12.

[0077] In this embodiment, the specific implementation process of selecting a suitable cluster to calculate the local Euclidean distance and local density is as follows:

[0078] Step 3.1: Condition A: If the number of clusters in the first clustering result is 1, then the k neighboring points are placed in the parameter calculation point cluster.

[0079] Step 3.2: Condition B: If the number of clusters in the first clustering result is 2, and the angle between the normal vector of the point with the largest elevation value among the k neighboring points and the z direction is less than the angle threshold T angle , then the farthest away from the detection point The clusters are put into the parameter calculation point clusters.

[0080] Step 3.3: If conditions A and B are not met, the nearest neighbor detection point The clusters are put into the parameter calculation point clusters.

[0081] Step 3.4: Calculate the local Euclidean distance of the points in the cluster according to the parameters. The specific calculation formula is:

[0082]

[0083] Among them, d Loc is the local Euclidean distance, k r Calculate the number of neighboring points in the cluster for the parameter, Point to be tested At the jth neighbor point in the reference phase point cloud, (x i ,y i ,z i ) is the point to be detected The three-dimensional coordinates, (x j ,y j ,z j ) is a point The three-dimensional coordinates of .

[0084] Step 3.5: Calculate the parameters to calculate the center of the cluster and search for k in the reference phase point cloud dThe local density is calculated based on the neighboring points. The specific calculation formula is:

[0085]

[0086] Among them, ρ Loc is the local density, k d is the number of neighboring points in the reference phase point cloud, r max Calculate the center of the cluster to k for the parameter d The maximum distance between neighboring points.

[0087] Step 4: Improve the local Euclidean distance according to the local density to obtain the density-adaptive local Euclidean distance. The specific calculation formula is:

[0088]

[0089] in, is the density-adaptive local Euclidean distance, d Loc is the local Euclidean distance, ρ Loc is the local density.

[0090] Step 5: All the points in the time phase to be detected whose density adaptive local Euclidean distance is greater than the distance threshold T d The points to be detected are extracted into the candidate change area, and the density adaptive local Euclidean distance is less than or equal to the distance threshold T d The points to be detected are considered unchanged and do not need to be extracted into the candidate change area.

[0091] In this embodiment, the distance threshold T d Set to 1.2m.

[0092] Step 6: Search for candidate change points q in the candidate change region i The k neighboring points on the xy plane in the phase point cloud to be detected are clustered by Euclidean method to obtain the second clustering result. The candidate change point q is clustered according to the second clustering result. i Perform processing to remove or not remove the candidate change area.

[0093] In this embodiment, the candidate change point q is clustered according to the second clustering result. i The specific implementation process of moving out or not moving out of the candidate change area is as follows:

[0094] Step 6.1: When the number of clusters in the second clustering result is 1, the candidate change point q i It does not belong to the area where the ground object modeling is poor, and the candidate change point q i Perform removal processing.

[0095] Step 6.2: When the number of clusters in the second clustering result is 2, 3 or 4, the candidate change point qi Belong to the area where the ground object modeling is not good, if the candidate change point q i The first top point of exists and has not changed, then the candidate change point q i Move out of the candidate change area; otherwise, the candidate change point q will not be i Processing; Among them, the candidate change point q i The first top point is the point where the angle θ between the normal vector of the k neighboring points and the z direction is less than the angle threshold T angle And the point with the largest elevation value.

[0096] Step 6.3: When the number of clusters in the second clustering result is greater than 4, the candidate change point q i Belong to the area where the ground object modeling is not good, if the candidate change point q i If the second top point of the candidate change point q changes, i Process; otherwise, the candidate change point q i Move out of the candidate change area; among them, the candidate change point q i The second top point of is the point with the largest elevation value among the k neighboring points.

[0097] Step 7: Cluster the points in the candidate change area processed in step 6, and calculate the area S corresponding to each cluster according to the average density of each cluster; if the area S is greater than the area threshold T area , then the cluster is considered to have truly changed; otherwise, the cluster is considered to be a false detection.

[0098] In this embodiment, the average density calculation formula of each cluster is:

[0099]

[0100] in, is the average density of the current cluster, k c is the number of points in the current cluster, (x j ,y j ,z j ) is the three-dimensional coordinate of the j-th point in the current cluster, (x n ,y n ,z n ) is the three-dimensional coordinate of the nearest point corresponding to the j-th point in the current cluster.

[0101] The calculation formula for area S is:

[0102]

[0103] Among them, S is the area of ​​the current cluster, is the average density of the current cluster, k cis the number of points in the current cluster. In this embodiment, the area threshold T area Set to 12 square meters.

[0104] by Figure 2 (b) is the reference time point cloud, Figure 2 (a) is the time phase point cloud to be detected (the first case). The detection results obtained by using the method of the present invention for three-dimensional change detection are as follows: Figure 3 (a) Figure 2 (a) is the reference phase point cloud, and (b) is the phase point cloud to be detected (the second case). The detection results obtained by using the method of the present invention for three-dimensional change detection are as follows: Figure 3 After the detection by the method of the present invention, each point to be detected in the detection result has a label (changed or unchanged), and then the detection result is compared with the Figure 4 The accuracy verification data of the manual marking (the label of each point to be detected is manually marked) is compared. If the two are consistent, the detection is considered correct, otherwise it is considered an error. The detection accuracy obtained is shown in Table 1:

[0105] Table 1 Detection accuracy of the three-dimensional change detection method of the present invention

[0106] category Precision Recall F1-score The first case 0.847 0.806 0.826 Second case 0.916 0.870 0.892

[0107] The traditional three-dimensional change detection method based on Euclidean distance segmentation is used for detection and compared with the manually labeled accuracy verification data. The detection accuracy is shown in Table 2:

[0108] Table 2 Detection accuracy of traditional three-dimensional change detection method based on Euclidean distance segmentation

[0109] category Precision Recall F1-score The first case 0.349 0.842 0.493 Second case 0.802 0.943 0.867

[0110] The evaluation metric uses the average F1 score to evaluate the detection performance. The F1 score is the harmonic mean of the precision and recall of each point. The formula is defined as follows:

[0111]

[0112] Where TP is the number of points for which the label detected by the present invention changes and the label manually marked for the point also changes; TN is the number of points for which the label detected by the present invention does not change and the label manually marked for the point also does not change; FP is the number of points for which the label detected by the present invention does not change but the label manually marked for the point changes; and FN is the number of points for which the label detected by the present invention changes but the label manually marked for the point does not change. As can be seen from Tables 1 and 2, the detection method of the present invention is far superior to traditional detection methods, with better noise immunity and higher detection accuracy.

[0113] The above disclosure is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with this technical field can easily think of changes or modifications within the technical scope disclosed in the present invention, and they should all be covered by the scope of protection of the present invention.

Claims

1. A method for detecting three-dimensional changes in ground objects, characterized in that: The following steps are involved: Step 1: Read the point clouds of two phases, perform precise alignment on them, and use one of the point clouds of the two phases after precise alignment as the reference phase point cloud, and the other as the phase point cloud to be detected; Step 2: Calculate the normal vectors of the time phase point cloud to be detected and the reference time phase point cloud; Step 3: Calculate each point to be detected in the time phase point cloud to be detected The angle θ between the normal vector and the z direction; If the angle θ is less than the angle threshold T angle , then search for the point to be detected k neighboring points on the xy plane in the reference phase point cloud, and clustering the k neighboring points to obtain a first clustering result, and selecting a suitable clustering cluster according to the first clustering result to calculate the local Euclidean distance and the local density; Otherwise search for the point to be detected k neighboring points in the three-dimensional space of the reference phase point cloud, and calculate the local Euclidean distance and local density; Step 4: improving the local Euclidean distance according to the local density to obtain a density-adaptive local Euclidean distance; Step 5: The density-adaptive local Euclidean distance in the phase point cloud to be detected is greater than the distance threshold T d The points to be detected are extracted into candidate change areas; Step 6: Search for candidate change points q in the candidate change region i k neighboring points on the xy plane in the phase point cloud to be detected, and clustering the k neighboring points to obtain a second clustering result, and the candidate change point q is clustered according to the second clustering result i Performing a process of moving out of or not moving out of the candidate change area; Step 7: Cluster the points in the candidate change area processed in step 6, and calculate the area S corresponding to each cluster according to the average density of each cluster; if the area S is greater than the area threshold T area , then the cluster is considered to have changed; otherwise, the cluster is considered to be a false detection; Preferably, in step 1, an ICP algorithm is used to accurately align the point clouds of the two phases.

2. The method for detecting three-dimensional changes of ground objects according to claim 1, characterized in that: In step 2, the specific implementation process of calculating the normal vectors of the time phase point cloud to be detected and the reference time phase point cloud is: Establishing a three-dimensional KD tree of the time phase point cloud to be detected and a three-dimensional KD tree of the reference time phase point cloud; The normal vectors of the to-be-detected time phase point cloud and the reference time phase point cloud are respectively calculated according to the three-dimensional KD tree.

3. The method for detecting three-dimensional changes of ground objects according to claim 1, characterized in that: In step 3, the calculation formula of the angle θ is: Among them, n x Point to be tested The value of the normal vector in the x direction, n y Point to be tested The value of the normal vector in the y direction.

4. The method for detecting three-dimensional changes of ground objects according to claim 1, characterized in that: In step 3, the specific implementation process of selecting a suitable cluster to calculate the local Euclidean distance and local density is as follows: Condition A: If the number of clusters in the first clustering result is 1, the k neighboring points are placed in the parameter calculation point cluster; Condition B: If the number of clusters in the first clustering result is 2, and the angle between the normal vector of the point with the largest elevation value among the k neighboring points and the z direction is less than the angle threshold T angle , then the point farthest from the detection point The clusters of are put into the parameter calculation point clusters; If conditions A and B are not met, the nearest neighbor detection point The clusters of are put into the parameter calculation point clusters; Calculate the local Euclidean distance of the points in the point cluster according to the parameters; Calculate the parameters to calculate the center of the point cluster and search for k in the reference phase point cloud d Calculate the local density of neighboring points; Preferably, the calculation formula of the local Euclidean distance is: Among them, d Loc is the local Euclidean distance, k r Calculate the number of neighboring points in the point cluster for the parameter, Point to be tested At the jth neighbor point in the reference phase point cloud, (x i ,y i ,z i ) is the point to be detected The three-dimensional coordinates (x j ,y j ,z j ) is a point The three-dimensional coordinates of Preferably, the calculation formula of the local density is: Among them, ρ Loc is the local density, k d is the number of neighboring points in the reference phase point cloud, r max Calculate the center of the point cluster for the parameter k d The maximum distance between neighboring points.

5. The method for detecting three-dimensional changes of ground objects according to claim 1, characterized in that: In step 4, the calculation formula of the density-adaptive local Euclidean distance is: in, is the density-adaptive local Euclidean distance, d Loc is the local Euclidean distance, ρ Loc is the local density.

6. The method for detecting three-dimensional changes of ground objects according to any one of claims 1 to 5, characterized in that: In step 6, the candidate change point q is clustered according to the second clustering result. i The specific implementation process of moving out or not moving out of the candidate change area is as follows: If the number of clusters in the second clustering result is 1, then the candidate change point q i It does not belong to the area where the terrain modeling is poor, and the candidate change point q i Carry out removal processing; If the number of clusters in the second clustering result is 2, 3 or 4, then the candidate change point q i It belongs to the area where the ground object modeling is not good. If the candidate change point q i If the first top point of exists and has not changed, then the candidate change point q i Remove candidate change regions; Otherwise, do not select the change point q i Processing; Among them, the candidate change point q i The first top point is the point where the angle θ between the normal vector of the k neighboring points and the z direction is less than the angle threshold T angle And the point with the largest elevation value; If the number of clusters in the second clustering result is greater than 4, the candidate change point q i It belongs to the area where the ground object modeling is not good. If the candidate change point q i If the second top point of the i Processing; otherwise, the candidate change point q i Move out of the candidate change area; where the candidate change point q i The second top point of is the point with the largest elevation value among the k neighboring points.

7. The method for detecting three-dimensional changes of ground objects according to any one of claims 1 to 5, characterized in that: In step 7, the average density calculation formula of each cluster is: in, is the average density of the current cluster, k c is the number of points in the current cluster, (x j ,y j ,z j ) is the three-dimensional coordinate of the jth point in the current cluster, (x n ,y n ,z n ) is the three-dimensional coordinate of the nearest point corresponding to the jth point in the current cluster.

8. The method for detecting three-dimensional changes of ground objects according to any one of claims 1 to 5, characterized in that: In step 7, the calculation formula of area S is: Among them, S is the area of ​​the current cluster, is the average density of the current cluster, k c is the number of points in the current cluster.

9. An electronic device comprising a memory and a processor, wherein the memory stores a computer program that can be run on the processor, wherein: When the processor runs the computer program, the processor executes the steps of the method for detecting three-dimensional changes of ground objects as described in any one of claims 1 to 8.

10. A computer-readable storage medium, wherein the computer-readable storage medium is a non-volatile storage medium or a non-transient storage medium, and a computer program is stored thereon, wherein: When the computer program is executed by a processor, the steps of the method for detecting three-dimensional changes of ground objects in any one of claims 1 to 8 are executed.

Images