Method and system for repairing missing areas of airborne lidar point cloud with existing DEM

Abstract

A method and system for repairing missing areas of airborne laser radar point cloud with reference to existing DEM, comprising: projecting the original point cloud to a two-dimensional plane, pre-processing the image using the closing operation in image morphology, clustering the points in the missing area to realize accurate positioning of the missing area, extracting the boundary points of the missing area using the eight-neighborhood detection method, filtering the boundary points to remove non-ground points, sampling the DEM according to the points in the missing area, and correcting the height under the constraint of the boundary points to obtain point cloud repair data. The image morphology is introduced into the positioning of the missing area, the height inverse distance weighted correction value of the existing DEM points preliminarily sampled in the missing area is constrained by the accurately extracted boundary points, and the repair of the missing area of the laser radar point cloud with high precision is realized. The present application better maintains the original terrain features of the missing area, and ensures the production of high-precision geographic information products with less cost.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of LiDAR data processing, and mainly aims at high-precision repair of airborne laser radar point cloud data with partial area missing, so as to meet the method of subsequent high-precision geographic product production requirements. BACKGROUND

[0002] With the continuous progress of LiDAR technology, more and more high-precision geographic products are applied. In the process of LiDAR data acquisition, due to the influence of ground object shielding, terrain undulation, water body, no-fly zone and the like, point cloud data cannot be obtained in local areas, and blank areas will appear in the point cloud, forming data missing. The missing area of the point cloud has no elevation information, and a high-precision ground model cannot be generated, which seriously affects the integrity and precision of the DEM product. The method of supplementing flight and reflight to fill the missing area of the point cloud needs to consume a large amount of manpower, financial resources and material resources, and the cost is too high. Therefore, a fast and economical internal processing method is needed to repair the missing area of the point cloud with high precision. The main methods for repairing the missing area at present are all from the point cloud itself, and mathematical interpolation and other methods are used to directly interpolate the three-dimensional point coordinates in the missing area. Commonly used are a point cloud hole repair method based on a boundary rejection algorithm proposed by Ming Ju [1] and others, which continuously pushes the boundary points into the hole to realize hole repair. This method is effective for small-range point cloud missing. A point cloud hole detection and repair method based on an ordered edge list method and a radial basis function is proposed by Hua Shuanggang [2] and others. This method uses a radial basis function to constrain on a regular surface to realize hole filling, but the repair effect of irregular hole areas needs to be improved. A hole filling algorithm based on local area is proposed by Changsheng Wang [4] and others, which uses a cubic spline interpolation method to construct a feature curve of the missing part. The consistency of the construction result of this method with the original terrain feature is slightly poor. A scattered point cloud model hole automatic repair method is proposed by Hao Linwen [5] and others, which is suitable for single hole conditions and does not consider the case where multiple holes exist. These methods cannot reflect the real terrain by only interpolating the elevation value of the point blank point according to the mathematical relationship without referring to the real spatial relationship. A hole repair method for laser radar point cloud data is proposed by Huang Zheng [5] and others. This method uses historical data as a reference for repair and considers the terrain feature, but the missing area detection method used by this method is greatly affected by the point cloud density, and the missed detection situation is easy to occur.

[0003] REFERENCES

[0004] [1]Ming Ju,Meili Wang.3D point cloud hole repair based on boundaryrejection method.2019,:105-108.

[0005] [2] Hua Shungang, Zhang Xiaoshuai. Point Cloud Hole Repair Combining Ordered Edge ListMethod and Radial Basis Function. 2020,09(01):56-61.

[0006] [3]Changsheng Wang,Ping Hu.A Hole-Filling Algorithm for TriangularMeshes in Engineering.2013,14(5):465-471.

[0007] [4] Hao Linwen. Research on hole repair technology in scattered point cloud model in reverse engineering [D]. Inner Mongolia University of Science and Technology, 2020.

[0008] [5] Huang Zheng, Wang Hongxing, Zhang Xin, Chen Yuquan, Wang Hainan. A vulnerability repair method for lidar point cloud data [P]. Jiangsu Province: CN114219917A, 2022-03-22. Summary of the Invention

[0009] This invention proposes a method and system for repairing missing areas of airborne lidar point clouds by referencing existing DEMs. It can achieve high-precision repair of missing areas of point clouds at extremely low cost. The repaired point clouds can better maintain the original terrain and have good practicality and universality.

[0010] 1. Locating the missing region:

[0011] Point cloud data is projected onto a two-dimensional plane. The two-dimensional image is preprocessed using the closing operation in image morphology to remove a large number of small noise gaps. The point-by-point clustering method is used to segment out the multiple missing regions. The boundary point cloud of each missing region is extracted using eight-neighbor detection. All boundary points are filtered to remove non-ground points.

[0012] 2. Missing region repair:

[0013] Calculate the elevation difference between the boundary points of each missing region and their corresponding coordinates on the existing DEM, and then calculate the average value T of the absolute value of the elevation difference. dBased on the points inside the missing regions, the existing DEM is sampled to obtain discrete points inside each missing region. The elevation correction value is then calculated using the inverse distance weighting method between each discrete point inside the missing region and its corresponding boundary points. The correction value is then converted to T... d After the constraints are applied and the discrete points within each missing region are corrected, the point cloud repair result is obtained.

[0014] The innovation of this invention lies in introducing image morphology into the localization of missing regions. By precisely extracting boundary points and constraining the elevation inverse distance weighted correction values ​​of existing DEM points initially sampled within the missing region, high-precision repair of missing regions in airborne LiDAR point clouds is achieved. This invention fully utilizes the potential of existing geographic data, better preserves the original terrain features of the missing region, ensures the production of high-precision geographic information products with minimal cost, and possesses excellent accuracy and universality.

[0015] This invention has the following characteristics:

[0016] 1. Accuracy: The accuracy of the results after repairing the missing areas of the simulated point cloud was evaluated. The accuracy met the national standards and could maintain a good fit with the point cloud existing in the original area and a certain correlation with the original terrain.

[0017] 2. Universality: This method can adapt to the repair of point cloud data of different densities and does not have high accuracy requirements for existing data containing elevation data. It can meet the needs of actual production.

[0018] 3. Economic efficiency: It can fully utilize the potential of existing geographic data and restore point cloud data at extremely low economic cost. Attached Figure Description

[0019] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings of the embodiments will be briefly described below.

[0020] Figure 1 A flowchart illustrating a method for repairing missing areas of airborne lidar point clouds using existing DEMs, as provided in an embodiment of the present invention.

[0021] Figure 2 Two sets of raw point cloud data are provided for one embodiment of the present invention.

[0022] Figure 3 The point cloud projection result is provided in one embodiment of the present invention.

[0023] Figure 4 The result of using an operator of size 9 to perform a closing operation and highlighting the missing region after accurate localization is provided in an embodiment of the present invention.

[0024] Figure 5The missing area repair result provided in one embodiment of the present invention (containing only ground points). Detailed Implementation

[0025] Figure 1 This paper demonstrates a method for repairing missing areas in an airborne lidar point cloud by referencing an existing DEM. For example... Figure 1 As shown, the method of the present invention includes the following 7 steps.

[0026] 1. Project the point cloud data onto a two-dimensional plane.

[0027] Based on the imaging principle of airborne point clouds, the average distance formula for point clouds is used:

[0028]

[0029] Where d is the average distance of the point cloud, and ρ is the point cloud density. The formula for the size of a 2D image is:

[0030]

[0031] Where H and W represent the height and width of the projected image, respectively; MAX and MIN are the maximum and minimum values ​​of the point cloud coordinates; X and Y represent the geodetic coordinates of the point cloud; and finally, the planar coordinate formula is used:

[0032]

[0033] To implement point cloud projection, 'a' represents projection using n times the average distance of the point cloud, which is usually taken as 1.

[0034] Figure 2 The projection results of the two sets of original point cloud data shown are available in [reference]. Figure 3 Besides the few clearly visible missing regions, the projected point cloud also contains a large amount of fine noise, which seriously affects the efficiency of subsequent clustering and segmentation.

[0035] 2. Use the closing operation in image morphology to preprocess the two-dimensional image to remove a large number of small noise gaps.

[0036] The semi-random scanning method results in uneven distribution of laser footprints, with some footprints having large gaps, leading to numerous small gaps in the projected point cloud image and affecting the detection of missing areas. In digital morphology, the process of dilation followed by erosion is called closing operation. Closing operation has the advantages of removing fine noise and smoothing boundaries. The dilation operation is defined as:

[0037]

[0038] The erosion operation is defined as follows:

[0039]

[0040] The closing operation is defined as:

[0041]

[0042] A represents the input image matrix, B represents the morphological operator, which is generally composed of an odd-numbered matrix of ≥3, and x represents the translation amount. It is the operator for expansion. It is the operator for erosion. Operators of different sizes will have a significant impact on the result of the closing operation. If the operator is too small, fine noise will not be sufficiently removed, while if the operator is too large, fine gaps will be easily enlarged.

[0043] Use operators of size 9 Figure 3 See the image for closing operations and highlighting of missing regions. Figure 4 ,contrast Figure 3 A large number of noise gaps were removed, and the missing area remained consistent with the original area. However, due to some operators with noise gaps larger than 9*9, a very small amount of noise was still difficult to remove. These noises could be removed by setting a threshold.

[0044] 3. The point-by-point clustering method was used to segment out the multiple missing regions.

[0045] 1) Determine a distance threshold T between categories and a gap missing threshold T2 based on prior knowledge.

[0046] 2) Determine a current point Pt as the first missing boundary point class.

[0047] 3) Calculate the distance Dis between the current point Pt and all the remaining boundary point clouds in sequence, and select the point Npt that is closest to Dis.

[0048] 4) If Dis is less than the threshold T, then add Npt to the class of the current point Pt and set Npt to Pt;

[0049] If Dis is greater than the threshold T, it means that the current missing region has been segmented and there are other missing regions. Add a new missing cluster and set Npt to Pt.

[0050] 5) Repeat steps 3-4 until all points are classified.

[0051] 6) Finally, based on the threshold T2, gaps that are misjudged as small missing regions and are difficult to remove due to gaps being larger than the closing operator are eliminated, and the missing region segmentation is completed.

[0052] 4. Eight-neighbor detection was used to extract the boundary point clouds of each missing region. All boundary points were then filtered to remove non-ground points. The results of this process can be found in [link to documentation]. Figure 4 .

[0053] 5. Calculate the elevation difference between the boundary points of each missing area and their corresponding coordinates on the existing DEM, and calculate the average value T of the absolute value of the elevation difference. d Calculation formula:

[0054] T d =∑|dH| / n

[0055] Calculate boundary constraint conditions T d T d The size of the value reflects the closeness between the Gpt and the DEM; the smaller the value, the better the closeness.

[0056] 6. Based on the points inside the missing regions, sample the existing DEM to obtain the discrete points inside each missing region.

[0057] 7. For each missing region, the elevation correction value is calculated using the inverse distance weighting method between discrete points within the region and their corresponding boundary points. The correction value is then weighted by T. d constraint.

[0058] There is an elevation difference dH between the ground point Gpt at the boundary of the missing region and the corresponding point on the DEM. This invention proposes the IDW-BC elevation correction method, which uses the inverse distance weighting method:

[0059]

[0060] When performing elevation correction, use the following formula:

[0061]

[0062] Constrain the inverse distance weighted correction value so that its magnitude is limited to [-T] d ,T d Within the range, if the correction value is outside the constraint value, then T with the same sign as ΔH is taken. d As the correction value, ΔH is the elevation correction value, dH is the elevation difference between the boundary point cloud and the DEM, D is the distance between the repaired point cloud and the boundary point cloud, and n is the number of boundary point clouds in the current missing region. The repaired point cloud result is obtained after correcting the discrete points within each missing region, as shown below. Figure 5 The results show that the boundary between the repaired point cloud and the original point cloud is smooth, and the terrain undulations are natural, resulting in a good repair outcome. This also proves the accuracy and universality of the invention.

[0063] In some embodiments, a system for repairing missing regions of lidar point clouds based on existing DEMs is also provided, comprising: a data acquisition module, a first processing module, a second processing module, a missing region segmentation module, a boundary point cloud extraction module for the missing region, a boundary constraint calculation module, an internal discrete point acquisition module, and a point cloud repair module.

[0064] The data acquisition module is configured to acquire point cloud data containing missing regions and existing DEM data containing elevation data for the corresponding regions. The first processing module is configured to project the point cloud data onto a two-dimensional plane to obtain a two-dimensional image. The second processing module is configured to preprocess the two-dimensional image using the closing operation in image morphology to remove small noise gaps. The missing region segmentation module is configured to segment the existing missing regions from the two-dimensional image after removing small noise gaps using a point-by-point clustering method. The missing region boundary point cloud extraction module is configured to extract the boundary point clouds of the missing regions using eight-neighbor detection and filter all boundary points to remove non-ground points. The boundary constraint calculation module is configured to calculate the elevation difference between the boundary points of the missing regions and their corresponding coordinates in the existing DEM data, and calculate the average value T of the absolute value of the elevation difference. d The internal discrete point acquisition module is configured to sample existing DEM data based on internal points within the missing regions to obtain internal discrete points for each missing region. The point cloud repair module is configured to calculate the elevation correction value ΔH using the inverse distance weighting method for each internal discrete point and its corresponding boundary point within the missing region. The elevation correction value ΔH is then used by T... d The point cloud restoration result is obtained by correcting the elevation correction values ​​for all discrete points within the missing region under constraints.

[0065] For more detailed implementation methods of each module of the system, please refer to steps 1-7 of the aforementioned method.

[0066] In some embodiments, a computer is also provided. The computer includes a processor and a memory. The memory is used to store non-transitory computer-readable instructions (e.g., one or more computer program modules). The processor is used to execute the non-transitory computer-readable instructions, which, when executed by the processor, can perform one or more steps of the method described above for repairing missing areas of an airborne lidar point cloud with reference to an existing DEM. The memory and the processor can be interconnected via a bus system and / or other forms of connection mechanisms.

[0067] For example, a processor can be a central processing unit (CPU), a graphics processing unit (GPU), or other form of processing unit with data processing and / or program execution capabilities. For instance, a CPU can be based on x86 or ARM architectures. A processor can be a general-purpose processor or a special-purpose processor, and it can control other components in a computer to perform desired functions.

[0068] For example, memory can include any combination of one or more computer program products, which can include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. Volatile memory can include, for example, random access memory (RAM) and / or cache memory. Non-volatile memory can include, for example, read-only memory (ROM), hard disk, erasable programmable read-only memory (EPROM), compact optical disc read-only memory (CD-ROM), USB storage, flash memory, etc. One or more computer program modules can be stored on the computer-readable storage medium, and the processor can run one or more computer program modules to implement various functions of the computer.

[0069] In some embodiments, a computer-readable storage medium is also provided for storing non-transitory computer-readable instructions that, when executed by a computer, can implement one or more steps of the method for repairing missing areas of an airborne lidar point cloud based on an existing DEM. That is, when the method and system for repairing missing areas of an airborne lidar point cloud based on an existing DEM provided in this application are implemented in software and sold or used as independent products, they can be stored in a computer-readable storage medium. For related descriptions of the storage medium, please refer to the corresponding description of memory in a computer above; further details will not be repeated here.

[0070] This invention introduces image morphology into the localization of missing regions. By precisely extracting boundary points, it constrains the elevation inverse distance weighted correction values ​​of existing DEM points initially sampled within the missing region, achieving high-precision repair of missing regions in airborne LiDAR point clouds. This invention fully utilizes the potential of existing geographic data, better preserves the original terrain features of the missing region, ensures the production of high-precision geographic information products with minimal cost, and possesses excellent accuracy and universality.

Claims

1. A method for repairing missing areas in a lidar point cloud by referencing an existing DEM, characterized in that, include: Acquire point cloud data containing missing regions and existing DEM data containing elevation data for the corresponding regions; Point cloud data is projected onto a two-dimensional plane to obtain a two-dimensional image; The closing operation in image morphology is used to preprocess the two-dimensional image to remove small noise gaps. The missing regions in a two-dimensional image with small noise gaps removed are segmented by point-by-point clustering, which includes: (1) determining a distance threshold T between categories and a gap missing threshold T2 based on prior knowledge; (2) Determine a current point Pt as the first missing boundary point class; (3) Calculate the distance Dis between the current point Pt and all the remaining boundary point clouds in turn, and select the point Npt with the closest Dis; (4) If the Dis is less than the threshold T, add Npt to the class of the current point Pt and set Npt to Pt; if the Dis is greater than the threshold T, add a new missing cluster and set Npt to Pt; (5) Repeat steps (3)-(4) until all points are classified; (6) Finally, remove gaps that are misjudged as small missing regions due to gaps being greater than the closing operator, based on the threshold T2, and complete the missing region segmentation; Eight-neighbor detection is used to extract the boundary point cloud of the missing region, and all boundary points are filtered to remove non-ground points; Calculate the elevation difference between the boundary points of the missing region and their corresponding coordinates in the existing DEM data, and then calculate the average value of the absolute elevation difference. ; The internal discrete points of each missing region are obtained by sampling the existing DEM data based on the internal points of the missing regions. The elevation correction values ​​were calculated by applying the inverse distance weighting method to the discrete points within the missing region and their corresponding boundary points. Elevation correction value quilt The point cloud restoration result is obtained after all discrete points within the missing region are corrected by elevation correction values, where the elevation correction values ​​are used to correct the point cloud. The size is limited to Within the range, if the elevation correction value Outside of the constraint value, take AND. Same symbol As the correction value, the constraint condition is as follows: 。 2. The method for repairing missing areas of lidar point clouds by referencing existing DEMs according to claim 1, characterized in that, A closed operation is achieved by using an odd matrix of ≥3 to perform dilation and erosion in image morphology on the projected two-dimensional image.

3. The method for repairing missing areas of lidar point clouds by referencing existing DEMs according to claim 1, characterized in that, The calculation formula is: Gpt represents the elevation difference between the ground point Gpt at the boundary of the missing region and the corresponding point in the existing DEM data, where n is the number of point clouds at the boundary of the missing region.

4. A system for repairing missing areas of a lidar point cloud based on an existing DEM, characterized in that, include: The data acquisition module is configured to acquire point cloud data containing missing regions and existing DEM data containing elevation data for the corresponding regions. The first processing module is configured to project point cloud data onto a two-dimensional plane to obtain a two-dimensional image. The second processing module is configured to preprocess the two-dimensional image using the closing operation in image morphology to remove small noise gaps. The missing region segmentation module is configured to segment existing missing regions from a two-dimensional image with small noise gaps removed using a point-by-point clustering method. It includes: (1) determining a distance threshold T between categories and a gap missing threshold T2 based on prior knowledge; (2) determining a current point Pt as the first missing boundary point class; (3) calculating the distance Dis between the current point Pt and all remaining boundary point clouds in sequence, and selecting the point Npt closest to Dis; (4) if Dis is less than the threshold T, adding Npt to the class of the current point Pt and setting Npt to Pt; if Dis is greater than the threshold T, adding a new missing cluster and setting Npt to Pt; (5) repeating steps (3)-(4) until all points are classified; (6) finally removing gaps that are misjudged as small missing regions due to gaps being greater than the closing operator, based on the threshold T2, thus completing the missing region segmentation. The missing region boundary point cloud extraction module is configured to extract the boundary point cloud of the missing region using eight-neighbor detection and filter all boundary points to remove non-ground points. The boundary constraint calculation module is configured to calculate the elevation difference between the boundary points of the missing region and their corresponding coordinates in the existing DEM data, and to calculate the average value of the absolute value of the elevation difference. ; An internal discrete point acquisition module is configured to sample existing DEM data based on points within the missing regions to obtain internal discrete points for each missing region; and The point cloud repair module is configured to calculate elevation correction values ​​for discrete points within the missing region and their corresponding boundary points using the inverse distance weighting method. Elevation correction value quilt The point cloud restoration result is obtained after all discrete points within the missing region are corrected by elevation correction values, where the elevation correction values ​​are used to correct the point cloud. The size is limited to Within the range, if the elevation correction value Outside of the constraint value, take AND. Same symbol As the correction value, the constraint condition is as follows: 。 5. The system for repairing missing areas of lidar point clouds by referencing existing DEMs according to claim 4, characterized in that, A closed operation is achieved by using an odd matrix of ≥3 to perform dilation and erosion in image morphology on the projected two-dimensional image.

6. The system for repairing missing areas of lidar point clouds by referencing existing DEMs according to claim 4, characterized in that, The calculation formula is: Gpt represents the elevation difference between the ground point Gpt at the boundary of the missing region and the corresponding point in the existing DEM data, where n is the number of point clouds at the boundary of the missing region.

Images