An airborne radar point cloud terrain product automatic production method and device

Abstract

The application discloses an airborne radar point cloud terrain product automatic production method and device, belongs to the field of image information processing, adopts a triangular net interpolation strategy, omits an interpolation range parameter, can produce a DSM / DEM terrain product without a cavity value, simultaneously supports a grid and point cloud data, overcomes input of multi-source data, processes each contour line in parallel, adopts a cubic Bezier curve to perform smooth fitting on the contour line, and efficiently and automatically produces a reasonable and smooth contour line product; in point cloud data processing, a gradual triangular net filtering algorithm is used as a basis to coarsely classify a terrain, a multi-level refined grid is established, and a complex terrain mixed scene is adapted; on the basis of a primary result, a point interval and a terrain growth angle are used as threshold values to perform secondary recovery on the terrain, and a situation of missing points of a ridge and a cliff terrain is made up; the algorithm rules a scene block, adopts a parallel processing strategy, and greatly improves filtering efficiency.

Description

Technical Field

[0001] This invention belongs to the field of image information processing, and in particular relates to an automated production method and apparatus for airborne radar point cloud terrain products. Background Technology

[0002] Airborne lidar (Light Detection and Ranging-LiDAR) is a novel active Earth observation technology. It overcomes the limitations of traditional point-by-point measurements, is unaffected by weather and lighting conditions, and can rapidly and cost-effectively acquire three-dimensional geographic information of large-scale surface scenes. Terrain products derived from this technology have wide applications in hydrology, meteorology, forestry surveys, power line inspection, geology, military, and national defense. Terrain product types mainly include: Digital Surface Models (DSM), Digital Elevation Models (DEM), and contour lines.

[0003] Since raw point cloud data does not contain classification information, point clouds of features such as buildings, vegetation, and bridges cannot accurately reflect the terrain within a scene. Therefore, automatically identifying ground and non-ground points from LiDAR point cloud data is a prerequisite for producing terrain products, i.e., point cloud filtering. Although many ground point filtering algorithms have been proposed, and many have achieved relatively ideal results in practical production applications, certain shortcomings still exist. First, existing algorithms cannot adaptively process data from complex scenes, requiring manual setting of different filtering parameters based on the scene's terrain features. For example, the Cloth Simulation Filter (CSF) (integrated in CloudCompare) requires manual specification of scene type parameters to complete filtering, resulting in poor filtering performance for scenes containing plains and mountains. The Progressive Morphological Filter algorithm completes filtering through morphological erosion and dilation operations, and the window size for erosion and dilation needs to be manually set based on the terrain and point density. The Progressive Triangulation Filter algorithm (integrated in Terrasolid) completes scene filtering by drawing a regular grid on the scene, taking the lowest point to construct the initial terrain, and gradually densifying it, but the selection of the initial grid size and the processing effect of terrain features such as gullies and ridges are unsatisfactory.

[0004] Digital surface model (DSM) and digital elevation model (DEM) are terrain products produced based on initial point clouds and ground points. Point cloud data is usually obtained by regular gridding. For ground point cloud data, the absence of non-ground points causes holes of different sizes in the ground point cloud. Common interpolation methods usually require manual specification of the interpolation range, which requires a certain amount of human intervention.

[0005] The production of contour map products is usually achieved through three-dimensional data acquisition, where elevation points are manually obtained at equal intervals and then interpolated to construct contour lines. Common automated extraction algorithms also have certain shortcomings, such as self-intersection of contour lines and jagged contour lines, which require post-processing and reduce production efficiency to some extent. Summary of the Invention

[0006] The present invention aims to solve the above problems and provides an automated production method and apparatus for airborne radar point cloud terrain products.

[0007] In a first aspect, the present invention provides an automated production method for airborne radar point cloud terrain products, including DSM / DEM grid terrain product production, DEM grid-based contour line production, and ground point-based contour line production.

[0008] The production of the DSM / DEM raster terrain product includes the following steps:

[0009] S11. First, process the point cloud data of the entire scene / ground.

[0010] S12. Construct a triangulation network for the processed point cloud data, determine whether to perform voxelization downsampling based on requirements, and perform triangulation network positioning pixel by pixel.

[0011] S13. Assign values ​​to the raster by interpolating the elevation value of the pixel position through the vertices of the triangular mesh; generate hole-free DSM / DEM data by interpolation through the establishment of a triangular mesh.

[0012] The production of contour lines based on DEM raster includes the following steps:

[0013] S21. First, set preset parameters, determine the elevation value of each contour line based on the input parameters, and generate binary classification mask data for the raster data using the specified elevation value. If the elevation value is higher than the specified value, mark it as 255; otherwise, mark it as 0. S22. Next, use linear interpolation to calculate contour points on the four sides of the raster and connect them to determine the type of contour line. To avoid contour lines passing through the corner points of the raster, increase the elevation value of all raster data by 0.00001 meters.

[0014] S23. Then, the contour point tracking method is used to track and process the contour points in both directions, starting from any contour line segment, to obtain a single contour line.

[0015] S24. Finally, smooth the contour lines to overcome the jagged effect.

[0016] The production of ground-point contour lines includes the following steps:

[0017] S31. First, set the preset parameters;

[0018] S32. Next, establish a triangulation network of ground points in the survey area and determine whether contour lines exist within the triangulation network. If so, use linear interpolation to determine contour points on the three sides of the triangulation network.

[0019] S33. Then, the contour point tracking method is used to track and process in both directions, starting from any contour line segment, to obtain a single contour line.

[0020] S34. Finally, the contour lines are smoothed to overcome the jagged effect.

[0021] Furthermore, in the automated production method of airborne radar point cloud terrain products described in this invention, step S11 specifically includes the following steps:

[0022] S111, Noise Removal: Noise removal is performed using statistical filtering. A Kdtree for point cloud data is established, and the standard deviation of the point cloud distance within the K-neighborhood is used. Noise data inevitably exists in the initial point cloud. It is usually assumed that the lowest point is the ground point, but in the actual data, there are many noise data that are much lower than the ground. These points usually exhibit discreteness and high variation.

[0023] S112, Scene coarse classification: Select the maximum building side length or the maximum non-ground point connected feature length in the scene to determine the grid size, select the lowest point in the grid as the grid value, solve the gradient in the xy direction, stretch the gradient image to 0-255, and use the Otsu's method to classify the scene into foreground and background values, and produce a binary map of the topographic relief of the survey area, where 255 is the foreground representing areas with large topographic relief and 0 is the background representing areas with small topographic relief.

[0024] S113, Progressive Triangulation Filtering: Initial ground point data is obtained through multi-level grid partitioning, and a triangulation network is established to restore the initial terrain. Newly added points are evaluated sequentially based on iteration angle, iteration distance, and mirror image judgment. If a point is a ground point, the triangulation network is updated; otherwise, it is skipped. The order of adding new points is based on the initial ground point, prioritizing those closest to it. Since points are generally closer to the initial ground point are more likely to be ground points, the order of adding new points directly affects the point cloud filtering effect.

[0025] S114. Secondary terrain restoration: Based on the results of the first filtering, ground points are used as the basis, and the second ground point judgment is performed by expanding outwards within a preset radius. The judgment conditions are point distance and growth angle, thus realizing the secondary restoration of terrain. Due to the uncertainty of the initial grid size and the defects of the triangulation judgment conditions, although mirror processing is introduced, ground points in locations such as gullies, steep slopes, and ridges often do not meet the judgment conditions of the triangulation. Therefore, secondary terrain restoration is performed. This step can restore the phenomenon of "flattened" ridges and defective boundaries of steep slopes.

[0026] S115, Misclassified Point Removal: Based on existing ground points, construct a DEM with a grid resolution of 2-3 times the average point interval. Select the lowest point in the point cloud within the grid as the grid value. Construct a buffer layer at a preset height (0.5 meters). Remove misclassified data above this buffer layer to obtain ground points and non-ground points, i.e., the processed point cloud data.

[0027] Furthermore, in the automated production method of airborne radar point cloud terrain products described in this invention, during the multi-level grid division process described in step S113, the mountainous area is refined into a secondary grid, with the grid size being half of the initial grid. The initial ground points in the area with large terrain undulations (i.e., the area marked as 255 in step S112) are supplemented, and error removal is performed simultaneously.

[0028] Furthermore, in the automated production method of airborne radar point cloud terrain products described in this invention, steps S113 and S114 are performed in a parallel processing manner; the size of a single processing unit is set, and the scene is divided into regular grid blocks; by adopting a parallel processing method, the efficiency of large scene point cloud traversal can be effectively improved, and the compatibility range of hardware memory can be increased.

[0029] Furthermore, in the automated production method of airborne radar point cloud terrain products described in this invention, virtual grid points are added to the segmented point cloud that forms a triangular network at the boundary of the grid segmentation.

[0030] Furthermore, in the automated production method of airborne radar point cloud terrain products described in this invention, the preset parameters include the initial elevation of contour lines, the final elevation, the basic contour interval, whether to produce calculation curves, and the type of auxiliary curve results.

[0031] Secondly, the present invention provides an automated production apparatus for airborne radar point cloud terrain products, comprising a processor and a memory electrically connected to each other; the memory is used to store computer programs; when the processor executes the aforementioned computer programs, it can realize the automated production method for airborne radar point cloud terrain products described in the first aspect.

[0032] Thirdly, the present invention provides a computer-readable storage medium storing a computer program; when the computer program is executed, it can realize the automated production method of airborne radar point cloud terrain products described in the first aspect.

[0033] The automated production method and apparatus for airborne radar point cloud terrain products described in this invention employs a triangular mesh interpolation strategy, omitting interpolation range parameters, and can produce DSM / DEM terrain products without holes or missing values. It simultaneously supports both raster and point cloud data, overcoming the challenges of multi-source data input, processing each contour line in parallel, and using cubic Bézier curves to smoothly fit the contour lines, efficiently and automatically producing reasonable and smooth contour line products. During point cloud data processing, a progressive triangular mesh filtering algorithm is used as the basis for coarse terrain classification, establishing multi-level refined grids to adapt to scenarios with complex mixed terrain. Based on the initial filtering results, a secondary terrain restoration is performed using point intervals and terrain growth angles as thresholds to compensate for missing points on ridges and steep slopes. The algorithm divides the scene into regular blocks and employs a parallel processing strategy, greatly improving filtering efficiency. Attached Figure Description

[0034] Figure 1 This is a schematic diagram of the point cloud data processing flow according to an embodiment of the present invention;

[0035] Figure 2 This is a schematic diagram of the automated production process of airborne radar point cloud terrain products according to an embodiment of the present invention;

[0036] Figure 3 This is a schematic diagram of a three-dimensional Kdtree as described in an embodiment of the present invention;

[0037] Figure 4 This is a schematic diagram of point cloud segmentation and initial ground point multi-level grid refinement as described in an embodiment of the present invention;

[0038] Figure 5 This is a supplementary schematic diagram of the multi-level grid refinement and encryption points and virtual points described in the embodiments of the present invention;

[0039] Figure 6 This is a schematic diagram of the mirroring process of the steep slope point according to an embodiment of the present invention;

[0040] Figure 7 This is a schematic diagram illustrating the definition of the terrain growth angle according to an embodiment of the present invention;

[0041] Figure 8 This is a schematic diagram of the iterative expansion process of secondary terrain restoration as described in an embodiment of the present invention;

[0042] Figure 9 This is a comparative schematic diagram of the filtering results described in the embodiments of the present invention;

[0043] Figure 10 This is a schematic diagram illustrating the possible existence of contour lines in the mask image according to an embodiment of the present invention;

[0044] Figure 11 This is a schematic diagram of cubic Bézier curve fitting according to an embodiment of the present invention;

[0045] Figure 12 This is a schematic diagram of the smoothing of the cubic Bézier curve of the contour lines according to an embodiment of the present invention;

[0046] Figure 13 This is a schematic diagram illustrating the presence of contour lines in the irregular triangular network described in an embodiment of the present invention.

[0047] Figure 14 This is a schematic diagram of bidirectional tracking of contour line segments in a triangular network according to an embodiment of the present invention;

[0048] Figure 15 This is a schematic diagram illustrating the terrain product results described in an embodiment of the present invention. Detailed Implementation

[0049] The automated production method and apparatus for airborne radar point cloud terrain products of the present invention will be described in detail below with reference to the accompanying drawings and embodiments.

[0050] Example 1

[0051] This disclosure provides an embodiment of a method for processing point cloud data of the entire scene / ground in an automated production method for airborne radar point cloud terrain products, such as... Figure 1 As shown, the specific process includes the following:

[0052] (1) Construct a point cloud Kdtree:

[0053] Point cloud data differs from raster data in that it is discrete data and cannot be used for neighborhood retrieval by index like raster data. Therefore, a Kd-tree needs to be built for the scene point cloud before processing to accelerate the neighborhood search speed. A Kd-tree is a data structure that partitions a k-dimensional data space. It is essentially a binary tree. Building a Kd-tree is a continuous partitioning process. Appropriate axes are selected as partitioning planes, and the depth is continuously increased until a Kd-tree for all point clouds is built. Point cloud retrieval is achieved by recursively searching and judging downwards from the root node.

[0054] (2) Point cloud noise removal:

[0055] In this embodiment of the disclosure, point cloud noise removal uses a distance standard deviation threshold between the point cloud and its neighboring points for removal, which is a global filtering algorithm. Figure 3 As shown, a point cloud Kdtree is established. The K-neighborhood of the current point is selected through the Kdtree. The distance between the current point and all neighboring points is calculated. The mean and standard deviation of the distances of all point clouds are calculated. The distance threshold is calculated using the following formula. Points greater than the threshold are considered noise and need to be removed.

[0056] Dis Max =Dis ave +nSigma*stddev

[0057] In the formula, Dis Max Dis is the distance threshold judgment condition. ave nSigma is the mean distance of all point cloud neighborhoods, nSigma is the standard distance deviation threshold that needs to be set manually, and stddev is the standard deviation of all point cloud neighborhood distances.

[0058] (3) Rough scene classification:

[0059] The length of the largest building or the largest connected feature in the scene is set, and the point cloud scene is divided into regular grids. The grid elevation value is selected based on the lowest elevation value falling into the grid. Considering that the point cloud may have hollow areas, interpolation processing is also required for the hollow areas. Gradient calculations are performed and merged in the x and y directions to obtain a comprehensive gradient map in the x and y directions. Gray-scale stretching is performed to stretch the gray values ​​to 0-255. The Otsu's method is used to perform binary classification on the gray-scale map. The foreground area is the area with obvious terrain undulations, and the background area is the area with gentle undulations. OTSU processes the gradient map instead of the depth map. Its purpose is to identify areas with obvious terrain changes (large gradients) rather than areas with high elevations.

[0060] The formula for calculating the gradient in the xy direction is:

[0061]

[0062] In the formula, gradf(x,y) is the gradient in the xy direction. These represent the gradient transformations in the x and y directions, respectively, which can be calculated by dividing the raster resolution by the elevation difference in the x and y directions.

[0063] The automatic determination method for the OTSU threshold can be achieved using the following formula.

[0064] D(T)=P B P O (μ B -μ O ),in In the formula, P B P O Let μ represent the probabilities of the background and the target, respectively. B μ O These represent the average gray values ​​for the background and target classes, respectively. T represents the threshold, L represents the gray level, N represents the number of pixels, and n represents the number of pixels at a certain gray level.

[0065] Figure 4 This is a diagram illustrating a rough classification of terrain. Figure 4 'a' represents the survey area scene, which mainly includes plains, hills, and gullies. Figure 4b represents the fused gradient map in the x and y directions. It can be observed that regions with gentle elevation changes have lower gray values, while regions with high gray values ​​show significant elevation changes. Figure 4 c represents the foreground and background regions in the Otsu's maximum inter-class difference method for automatic threshold classification, with areas exhibiting drastic slope changes being marked.

[0066] (4) Point cloud segmentation and initial ground point multi-level grid refinement:

[0067] Based on the obtained OTSU binary classification results, the foreground (i.e., grayscale of 255) is refined and densified by grid. The size of the refined grid is half that of the initial grid, and the grid value is still the lowest value of the point cloud within the grid. Newly added points also need to be removed as coarse errors. The point cloud is divided into regular grid blocks, and each block is a parallel processing unit. Considering the problem of building the grid at the boundary of the survey area and the block point cloud, virtual point cloud supplementation is required. The elevation of the supplemented points is selected from the elevation of the nearest point. Figure 5 This is a schematic diagram of virtual points added to the multi-level grid for refining densification points and foreground areas. Compared to regular grid points, adding densification points in areas with drastic terrain undulations can improve the accuracy of terrain description. Due to the small size of this scene, the segmentation effect is not shown at this time.

[0068] (5) Asymptotic triangular mesh filtering:

[0069] An irregular triangulated network (TIN) is built using initial ground points. The distances between the point cloud within the scene and the initial ground points are calculated and sorted. Starting with the minimum distance, each point is checked point by point to see if it meets the encryption conditions, which are the iteration angle and iteration distance. First, it's determined which triangle the current point falls within in the XOY plane of the projected triangulation. Then, the angles between the current point and the three vertices of the triangulation, as well as the distance from the point to the triangulation plane, are calculated and compared with the iteration angle and iteration distance thresholds. If the distance is less than the threshold, the point is considered a ground point, and the ground triangulation is updated accordingly. To improve the efficiency of network construction and triangulation, a minimum side length for the triangulation can be set to control the updates.

[0070] The progressive triangulation filtering algorithm incorporates a mirroring process for handling steep slopes. If the angle between the triangulation and the XOY plane reaches a preset angle, it is considered a steep slope point, i.e., P in the figure below. potential Find the nearest ground point P, and calculate the mirror point P. mirror For the mirror point, calculate the iterative angle and iterative distance. If the condition is met, then P... potential This method recovers some ground point data to a certain extent, but the processing effect is not satisfactory in actual engineering applications.

[0071] (6) Secondary terrain restoration

[0072] In actual production, it was found that many ground points, such as steep slopes, gullies, and ridges, could not be recovered using the above method. The reason for this was that these points did not meet the conditions for triangulation. Therefore, a second recovery of the existing terrain was performed based on the first filtering. The terrain growth angle is defined as the angle between the current point and the fitted plane within the R-neighborhood of the nearest ground point. Figure 7 The term θ represents the topographic angle. The specific method involves iteratively expanding existing ground points within a preset radius. For each expanded point, the terrain growth angle needs to be calculated; if it is less than a threshold, it is included as a ground point. Figure 8 The yellow area shown represents the result of the initial filtering of the progressive triangulation network, while the blue, red, green, and pink areas represent the results of iterative expansion during the secondary terrain restoration process. The point cloud missing in the steep slope boundary area has been restored.

[0073] (7) Removal of misclassified points

[0074] The method described in this embodiment inevitably generates some misclassified points after processing the scene point cloud, especially during the second terrain restoration process, where some non-ground points also meet the terrain growth conditions. Therefore, it is necessary to remove these misclassified points. A DEM is constructed based on existing ground points, with the grid value selected as the lowest point cloud value within the grid. A buffer layer with a preset height (0.5 meters) is then built on the DEM. Points above this buffer layer need to be classified as non-ground points.

[0075] (8) Ground point filtering ends

[0076] The process described above can adaptively and efficiently complete the filtering of ground points and non-ground points in complex scenarios, and this data provides data support for the production of subsequent terrain products.

[0077] For the data in this test area, in this embodiment of the disclosure, the filtering algorithms built into PCM and CloudCompare software are compared with the filtering results of the present invention. Figure 9 As shown, 9a is a schematic diagram of the triangular mesh filtering result; Figure 9 b is a schematic diagram of the morphological filtering results; Figure 9 c is a schematic diagram of the fabric simulation filtering results; Figure 9 d is a schematic diagram of the processing result of the method described in this embodiment; red represents non-ground points and blue represents ground points. It can be clearly seen that the filtering result of the method described in this embodiment is the best, and it can achieve relatively satisfactory processing results in complex terrain conditions, resulting in processed point cloud data.

[0078] Example 2

[0079] Based on the above embodiment one, this disclosure presents an automated production method for airborne radar point cloud terrain products, such as... Figure 2The process shown includes the following steps:

[0080] (1) DSM and DEM production

[0081] First, determine the resolution of the output image. If not specified, the resolution is output based on the average point interval. The average point interval is determined by calculating the average distance between each point and its nearest neighbor.

[0082] Based on the scene point cloud data, the bounding box range is calculated, and the affine transformation matrix of the raster image is determined accordingly. The position of the raster product in the upper left plane of the bounding box is determined. The image width and height are calculated by the length and width of the measurement area and the image resolution. All initial image raster values ​​are set to -9999 invalid values.

[0083] Based on the resolution, determine whether voxel downsampling of point cloud data is necessary. That is, draw regular grids in the x, y, and z directions, with the downsampling interval being the image resolution. Select the point within a voxel that is closest to the voxel center as the voxel sampling point. This method can reduce the number of triangulations to be constructed and improve the algorithm's running efficiency without affecting the raster accuracy.

[0084] An irregular triangular network (TIN) is constructed from the downsampled point cloud data. Pixels are traversed row by row and column by column, converting pixel coordinates to geographic coordinates to determine which TIN the pixel belongs to. Based on the spatial coordinates of the three vertices of the TIN, a plane equation is fitted, and the Z-value coordinates are obtained from the plane coordinates and assigned to the raster image. The methods for converting pixel coordinates to geographic coordinates, determining whether a point is inside a triangle (arbitrary polygon), and calculating the spatial plane equation are as follows:

[0085] XGeo=GeoTransform[0]+GeoTransform[1]*X pixel +Y pixel *GeoTransform[2]

[0086] YGeo=GeoTransform[3]+GeoTransform[4]*X pixel +Y pixel *GeoTransform[5]

[0087] In the formula, the geographic coordinates XGeo and YGEO are respectively represented by X p ixe,lY pixel Pixel coordinates are calculated using the affine transformation matrix GeoTransform.

[0088] t=(Px-Bx)*(Ay-By)-(Ax-Bx)*(Py-By)

[0089] In the formula, A and B represent the previous and next points of the polygon, respectively, x and y represent the planar coordinates, P represents the current point, and the formula iterates through all points. If t < 0 and the number of points is odd, then the point is inside any polygon; otherwise, it is not.

[0090] Ax + By + Cz + D = 0 (C ≠ 0)

[0091]

[0092] make:

[0093]

[0094] Solve for x using least squares: x=(A T A) -1 A T b

[0095] In the formula, (x i ,y i ,z i ) represents the coordinates of a three-dimensional point involved in calculating the equation of the spatial plane, where A, B, C, and D are the coefficients of the plane equation.

[0096] (2) Contour production

[0097] To ensure compatibility with contour line production from different data sources, this embodiment provides contour line production strategies based on DEM raster and contour line production strategies based on ground point cloud data. This approach can meet the diverse input of data and improve the adaptability of the algorithm.

[0098] The contour line production process based on DEM raster includes:

[0099] First, confirm the starting and ending elevations of the contour lines, the basic contour interval, and whether they are production curves and auxiliary curves, etc. Calculate the number of contour lines that need to be produced in the scene and their corresponding elevations, and process them in parallel with each contour line at each elevation as the smallest unit.

[0100] For each contour line, a mask image of the same size as the raster image is created and initialized to 0. All DEM rasters are traversed. If the pixel value of a raster is greater than the elevation value of the contour line, the pixel is assigned a value of 255, and a mask image of the corresponding contour line is generated.

[0101] Figure 10This represents all possible scenarios in the masked image. Black values ​​correspond to a mask image value of 0, white values ​​correspond to a mask image value of 255, and red lines represent contour lines. Pixel-by-pixel analysis is performed. If a pixel is the first or last in the image below, it is not processed. If it belongs to one of the remaining 14 grid types, a corresponding contour line segment is generated for that grid. The endpoints of the contour line segment must intersect the four sides of the grid at two / four points (one / two contour lines). The endpoint positions are determined using linear interpolation. Considering that endpoint values ​​may be located at grid corners, each grid value is incremented by 0.00001 meters.

[0102] In this embodiment, contour line tracking includes: setting all pixels of the mask to an unmarked state, and marking the grid with two contour lines; starting from any contour line, recording the positions of its two endpoints, searching for the next contour line from one endpoint, marking the traversed grids until all contour lines are found, and then searching in reverse from the other endpoint, repeating the same method until all segments of the contour line are found, and connecting them sequentially in reverse at the endpoints; if the two endpoints of a contour line are close together, the contour line is closed; repeating this process can obtain all contour lines for the current elevation value; if a contour line is too short, a distance threshold can be used to delete the contour line to avoid the occurrence of local closure of multiple contour lines due to terrain anomalies.

[0103] In this embodiment of the invention, contour smoothing processing includes the following: Contour lines produced based on a raster are actually formed by sequentially connecting multiple contour segments. The resulting contour lines will appear "jagged". To remove this phenomenon, the present invention uses cubic Bézier curves to smooth the contour lines, and the produced contour lines are more in line with the actual survey area scene. It is worth mentioning that the cubic Bézier curve is formed by reconstructing and adding contour points with a preset step size (such as 0.1) using four contour points as control points and connecting them sequentially. This refines the contour point data to a certain extent, and can also fit more refined contour line data for raster images with higher resolution.

[0104] Bézier curves are actually controlled by n control points to define the shape of the curve. They are typically used to draw vector curves. Bézier curves pass through the first and last two points, but not through the midpoint. By setting an appropriate step size, ... Figure 11 Taking the four contour points P0, P1, P2, and P3 as an example, the route traversed by the curve is as follows.

[0105] The cubic Bézier curve is calculated as follows:

[0106]

[0107] In the formula, B(x,t), B(y,t), and B(z,t) represent the three-dimensional spatial coordinates of the Bézier curve points calculated at step size t using four contour control points P0, P1, P2, and P3.

[0108] Figure 12 This illustrates the practical application of cubic Bézier curves in smoothing contour lines. Black represents the contour lines before smoothing, while red represents the smoothed result. It is easy to see that the smoothed contour lines overcome the "jagged" phenomenon, are more aesthetically pleasing, and better reflect the actual terrain changes.

[0109] In this embodiment of the disclosure, the contour lines are vectorized, and the vector type is polyline. The attribute fields are Shape (vector type), ID (contour line ID), elev (elevation field) and Type (contour line type: primary curve, auxiliary curve, calculation curve and intermediate curve).

[0110] The process of generating contour lines based on ground point clouds includes:

[0111] First, confirm the starting and ending elevations of the contour lines, the basic contour interval, and whether they are production curves and auxiliary curves, etc. Calculate the number of contour lines that need to be produced in the scene and their corresponding elevations, and process them in parallel with each contour line at each elevation as the smallest unit.

[0112] An irregular triangular mesh is created from ground point cloud data. Unlike raster imagery, in a triangular mesh, contour lines only exist in certain locations. Figure 13 In the three cases shown, the endpoint positions are obtained by linear interpolation. Considering that the endpoint values ​​may be located at the vertices of the triangulation network, the operation of increasing the elevation of each triangulation network vertex by 0.00001 meters can effectively avoid this. Contour segments with corresponding elevations are generated in the triangulation network respectively.

[0113] The following formula can be used to calculate and determine whether contour segments exist in a triangular network:

[0114]

[0115] If n1+n2+n3=2, it means that contour segments exist in the triangular network.

[0116] In the formula, P1, P2, and P3 are the coordinates of the three vertices of the triangulation network, .z represents the z-coordinate of the corresponding point, and P represents the elevation of the contour line.

[0117] Contour segment tracing and contour smoothing have been described in detail in the aforementioned description of contour generation based on DEM raster, and will not be repeated here. The methods are the same as above. Figure 14This diagram illustrates bidirectional tracing of contour lines in an irregular triangular network. Red represents the initial contour line, yellow represents the contour line traced from the left endpoint, and blue represents the contour line traced from the right endpoint. Therefore, it is necessary to preserve the relationship between the triangular network and adjacent triangular networks during the construction of the network. The contour line vectors are written out, and the content and method are the same as above.

[0118] Based on the above process in this embodiment, various terrain products are efficiently and automatically produced for this survey area. Some terrain product results are as follows: Figure 15 As shown, where Figure 15 a is a schematic diagram of the DSM topographic product for the survey area. Figure 15 b is a schematic diagram of the DEM topographic product of the survey area. Figure 15 c is a schematic diagram of the topographic contour products of the survey area. Figure 15 d is a partial contour map of the mountain and cliff; the DSM / DEM resolution is 3 meters for the average point interval of the point cloud, and the contour line products have a basic contour interval of 10 meters. Since the contour lines produced based on the raster and the point cloud are almost overlapping, only one product is shown.

[0119] In terms of algorithm efficiency, with approximately 96,000+ point clouds in the survey area and hardware conditions including an 8-core CPU and 16GB of RAM, point cloud filtering took 2 seconds, DSM / DEM generation took 1 second, and contour line generation took 1 second, indicating relatively high algorithm efficiency. From a product quality perspective, the DSM / DEM product can realistically reflect the terrain features and landforms in the scene; the contour line product closely adheres to the point cloud surface, accurately reflecting the topographic undulations of the survey area, and performs well in plains, mountains, and cliffs.

[0120] Example 3

[0121] This disclosure presents an automated production device for airborne radar point cloud terrain products, comprising a processor and a memory electrically connected to each other; the memory is used to store computer programs; when the processor executes the aforementioned computer programs, it can implement the automated production method for airborne radar point cloud terrain products described in Embodiments 1 and 2 above. The specific production method is the same as that in Embodiments 1 and 2 above, and will not be repeated here.

[0122] Example 4

[0123] This disclosure provides a computer-readable storage medium storing a computer program. When the computer program is executed, it can implement the automated production method of airborne radar point cloud terrain products described in Embodiments 1 and 2 above. The specific production method is the same as that in Embodiments 1 and 2 above, and will not be repeated here.

[0124] The computer described in this application embodiment can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable devices. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. The computer-readable storage medium can be any usable medium that a computer can read, or a data storage device such as a server or data center that integrates one or more usable media. The usable medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., digital versatile optical disc (DVD)), or a semiconductor medium (e.g., solid-state drive (SSD)). The software formed by the computer's stored code can reside in random access memory, flash memory, read-only memory, programmable read-only memory, electrically erasable programmable memory, registers, or other storage media that are mature in the art.

[0125] In the various embodiments of this application, the functional modules can be integrated into one processing unit or module, or each module can exist physically separately, or two or more modules can be integrated into one unit or module. In the above embodiments, they can be implemented entirely or partially through software, hardware, firmware, or any combination thereof. When implemented using software, they can be implemented entirely or partially in the form of a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated.

[0126] The automated production method and apparatus for airborne radar point cloud terrain products described in this invention can adaptively, efficiently, and automatically produce airborne LiDAR point cloud ground point filtering, DSM / DEM, and contour terrain products in complex scenarios. The algorithm exhibits strong data compatibility, adaptability to complex terrain, and minimal parameter adaptability. Extensive data experiments and engineering production tests have demonstrated that this method and apparatus perform well.

Claims

1. An automated production method for airborne radar point cloud terrain products, characterized in that: This includes the production of DSM / DEM raster terrain products, the production of DEM raster contour lines, and the production of ground point contour lines. The production of the DSM / DEM raster terrain product includes the following steps: S11. First, process the point cloud data of the entire scene / ground. Step S11 specifically includes the following steps: S111, Noise Removal: Noise removal is performed using statistical filtering. A point cloud data Kdtree is established, and the standard deviation of the point cloud distance within the K-neighborhood is used. S112, Scene coarse classification: Select the maximum building side length or the maximum non-ground point connected feature length in the scene to determine the grid size, select the lowest point in the grid as the grid value, solve the gradient in the xy direction, stretch the gradient image to 0-255, and use the maximum inter-class difference method to classify the scene into foreground and background values, and produce a binary map of the topographic relief of the survey area, where 255 is the foreground representing areas with large topographic relief and 0 is the background representing areas with small topographic relief. S113, Progressive triangulation filtering: The initial ground point data is obtained by multi-level grid division and a triangulation network is established to restore the initial terrain. The newly added points are judged sequentially based on the iteration angle, iteration distance, and mirror judgment. If it is a ground point, update the triangulation; otherwise, skip it. The order of adding new points is based on the initial ground point, with those closer to it being prioritized. S114. Secondary terrain restoration: Based on the results of the first filtering, ground points are used as the basis, and the ground points are judged in a second round by expanding outwards within a preset radius. The judgment conditions are the point distance and the growth angle, so as to realize the secondary terrain restoration. S115, Removal of misclassified points; Based on existing ground points, a DEM is constructed with a grid resolution of 2-3 times the average point interval. The grid value is selected as the lowest point in the point cloud within the grid. A buffer layer is constructed with a preset height. Misclassified data above the buffer layer are removed to obtain ground points and non-ground points, i.e., the processed point cloud data. S12. Construct a triangulation network for the processed point cloud data, determine whether to perform voxelization downsampling based on requirements, and perform triangulation network positioning pixel by pixel. S13. Assign values ​​to the raster by interpolating the elevation value of the pixel position through the vertices of the triangular mesh; generate hole-free DSM / DEM data by interpolation through the establishment of a triangular mesh. The production of contour lines based on DEM raster includes the following steps: S21. First, set the preset parameters, determine the elevation value of each contour line based on the input parameters, and generate binary classification mask data for the raster data based on the specified elevation value. If the elevation value is higher than the specified value, mark it as 255; otherwise, mark it as 0. S22. Next, use linear interpolation to calculate contour points on the four sides of the grid and connect them to determine the type of contour lines; increase the elevation value of all grids by 0.00001 meters. S23. Then, the contour point tracking method is used to track and process the contour points in both directions, starting from any contour line segment, to obtain a single contour line. S24. Finally, smooth the contour lines. The production of ground-point contour lines includes the following steps: S31. First, set the preset parameters; S32. Next, establish a triangulation network of ground points in the survey area and determine whether contour lines exist within the triangulation network. If so, use linear interpolation to determine contour points on the three sides of the triangulation network. S33. Then, the contour point tracking method is used to track and process in both directions, starting from any contour line segment, to obtain a single contour line. S34. Finally, smooth the contour lines.

2. The automated production method for airborne radar point cloud terrain products according to claim 1, characterized in that: In step S113, during the multi-level grid division process, the mountainous area is refined into a secondary grid, with the grid size being half that of the initial grid. Initial ground points in areas with large terrain undulations are supplemented, and error removal is performed simultaneously.

3. The automated production method for airborne radar point cloud terrain products according to claim 2, characterized in that: Steps S113 and S114 are performed in parallel; the size of a single processing unit is set, and the scene is divided into regular grid blocks.

4. The automated production method for airborne radar point cloud terrain products according to claim 3, characterized in that: The grid segmentation adds virtual grid points to the point cloud of the triangular mesh at the boundary.

5. The automated production method for airborne radar point cloud terrain products according to claim 4, characterized in that: The preset parameters include the initial elevation of the contour line, the final elevation, the basic contour interval, whether to produce a curve, and the type of auxiliary curve result.

6. An automated production device for airborne radar point cloud terrain products, comprising a processor and a memory electrically connected in phase; the memory being used to store computer programs; characterized in that: When the processor executes the aforementioned computer program, it can realize the automated production method of airborne radar point cloud terrain products as described in any one of claims 1-5.

7. A computer-readable storage medium, characterized in that: The storage medium stores a computer program; when the computer program is executed, it can realize the automated production method of airborne radar point cloud terrain products as described in any one of claims 1-5.

Images