Digital elevation model generation method based on morphological filtering iteration

Abstract

The invention discloses a morphological filtering iteration-based digital elevation model generation method. The method comprises the following steps of: 1) partitioning an image; 2) morphological grayscale filtering; (3) iterative binarization is carried out; 4) marking connected domains; 5) threshold filtering; 6) surface interpolation; and 7) merging processing. According to the method, morphological gray filtering and iterative threshold segmentation technologies are fused, a block optimization strategy and a dynamic connected domain screening mechanism are combined, and the problems that a traditional algorithm is sensitive in window size, poor in terrain adaptability, multiple in set parameters and the like are effectively solved. And the automatic processing efficiency of generating a digital elevation model (DEM) by a digital surface model (DSM) under a complex terrain and the terrain restoration precision are obviously improved.

Description

Technical Field

[0001] The present invention relates to digital elevation model (DEM) generation technology and belongs to the field of geographic information processing technology. Specifically, it is a method for generating a digital elevation model based on morphological filtering iteration, which can efficiently and accurately filter out ground objects of the digital surface model and generate a DEM. Background Art

[0002] The digital elevation model (DEM) is an important data model reflecting the surface terrain undulation and is widely used in fields such as geological disaster prediction, water conservancy project construction, and environmental protection. Currently, the generation of DEM usually relies on digital surface model (DSM) data. However, the DSM contains ground objects such as buildings and vegetation, and these ground objects need to be removed to obtain an accurate DEM.

[0003] In the field of digital terrain model processing, for the problem of converting a digital surface model (DSM) into a digital elevation model (DEM), the traditional method is manual editing. Ground object points are screened and removed through empirical discrimination. This method requires a large amount of labor and time, and requires experienced workers to ensure accuracy.

[0004] With the progress of computer image automatic processing, many DSM processing algorithms have been applied, and a slope-based filtering algorithm has emerged. This algorithm believes that the slope value between a ground point and a non-ground point is greater than the slope value between ground points, and thus screens out the ground object points in the DSM. The advantage is high speed and efficiency. However, a single fixed slope threshold makes it perform poorly in mountainous areas with complex terrain.

[0005] There is also morphological filtering. By applying morphological opening operation, the height change of corresponding points before and after the opening operation is compared. Points exceeding the threshold are considered ground object points. However, the algorithm depends on the size of the sliding window. If it is too small, it cannot filter out large-area ground objects. If it is too large, it will filter out undulating terrain.

[0006] Subsequently, progressive morphological filtering was proposed to address the shortcomings of the morphological algorithm. Sliding windows of different sizes are used for morphological operations, and finally the results of different sizes are weighted and calculated to obtain the final filtering result. This better solves the problem of difficult window size determination, but there will still be a phenomenon that the top of the area with a large slope is over-filtered.

[0007] The progressive irregular triangulation network filtering algorithm selects initial ground seed points to construct a triangulation network. Through conditional judgment, unclassified points are added to the triangulation network, and the iteration continues until no new points are added to the triangulation network. The points constituting the triangulation network are used to generate a DEM to achieve the effect of removing ground objects. This algorithm has high filtering accuracy, but depends on the selection of initial seed points and has low iteration efficiency.

[0008] In recent years, cloth simulation filtering has been proposed based on the principle of physical simulation. This algorithm flips and inverts all points. Assuming that a rigid cloth covers the point plane, the surface formed after the cloth sinks due to gravity is the required DEM plane. The cloth filtering algorithm has the advantages of few parameters and good filtering effect. However, different cloth parameters need to be set for different terrain conditions, and it performs poorly in areas with complex ground object types. Summary of the Invention

[0009] To solve the defects of the above existing algorithms in the process of generating DEM, the present invention provides an iterative algorithm based on morphological filtering, which can effectively remove ground objects from DSM and generate a high-precision digital elevation model (DEM). This method obtains the terrain surface through a morphological algorithm, iteratively increases the height value of the terrain surface to cut and obtain a sequence of binary images of the ground surface, determines the ground object range based on the elevation difference between the edges extracted from the binary image sequence and the pixels of the original image, and then filters out the ground objects and performs interpolation, thereby realizing the filtering of ground objects in DSM and the reconstruction of the DEM ground surface.

[0010] The object of the present invention is achieved as follows:

[0011] A method for generating a digital elevation model based on morphological filtering iteration, comprising the following steps:

[0012] 1) Image block division:

[0013] The input digital surface model DSM image is divided into image blocks according to a certain size (each image block is processed separately in the subsequent steps, and finally the processing results are merged). As an optimization, an overlapping area is set between blocks, and the pixel values in the overlapping area use the graph cut algorithm to find the splicing gaps between blocks.

[0014] Block processing is a common method. For large-scale remote sensing images with tens of millions or hundreds of millions of pixels, by loading and processing block by block, it can not only improve memory efficiency and avoid the problem of insufficient device memory, but also reduce the computational complexity and the overall processing difficulty.

[0015] When setting the size of the blocks for block processing in the invention, the size of each image block is set to be larger than the maximum size of the ground objects in the area, that is, the number of horizontal and vertical pixels of the image block is greater than the maximum number of horizontal and vertical pixels of the ground object contour. At the same time, in order to make the transition between image blocks smooth, a certain overlapping area is set between blocks, and the pixel values in the overlapping area use the graph cut algorithm to find the reasonable splicing gaps between blocks. The graph cut algorithm is used to optimize the block-to-block splicing to ensure terrain continuity.

[0016] Graph Cut is an image segmentation algorithm based on the minimum energy path. It regards the image as a graph and finds the optimal stitching boundary through Min-Cut to achieve smooth fusion:

[0017] M(s, t, A, B) = ∥A(s) - B(s)∥ + ∥A(t) - B(t)∥

[0018] s and t are two adjacent pixels in the overlapping region. A(s) and B(s) represent the pixel values of pixel s in image A and image B, and A(t) and B(t) represent the pixel values of pixel t in image A and image B. The purpose of the algorithm is to find a cut seam to minimize ∑ for all s,t on t he cut M(s, t, A, B), and the pixels on both sides of the cut seam take values from image A and image B respectively.

[0019] 2) Morphological gray-scale filtering:

[0020] Perform the opening operation and erosion operation of morphological operations on the image block in step 1). Use the structural element to perform the opening operation on the DSM image to obtain a rough terrain surface, and then use the structural element erosion operation to obtain the DSM with the edges of the objects eroded. Use the rough terrain surface as the threshold base map for iterative binarization, and use the DSM with the edges of the objects eroded as the base map for height difference calculation.

[0021] Morphological gray-scale filtering is an extension of mathematical morphology and is mainly used to process gray-scale images. The present invention uses the opening operation and erosion operation of morphological gray-scale filtering:

[0022]

[0023] O(x, y) = Dilation(Erosion(I(x, y)))

[0024] I(x, y) is the original image, E(x, y) is the erosion operation, O(x, y) is the opening operation, and S is the structural element.

[0025] For further optimization, perform an opening operation on the DSM image using a suitable structural element to obtain a rough terrain surface. In the opening operation, the size of the structural element needs to be larger than the radius of the largest feature in the area; otherwise, there will be feature protrusions in the rough terrain, affecting subsequent iterative binarization. As large a structural element as possible should be used. Then, perform an erosion operation to obtain the DSM with the edges of the eroded features. In this step, the structural element should not be set too large. Setting the size of the structural element (filter) to be less than 7×7 and smaller sizes is sufficient. Use the two results as the threshold base map for iterative binarization and the base map for height difference calculation, respectively.

[0026] 3) Iterative binarization:

[0027] To screen out features, a binary height threshold that can separate features from the ground needs to be found, and the image is binarized using this threshold to segment ground points and non-ground points.

[0028] The present invention uses the method of iterative binarization. For the image block obtained in step 1), the method of iterative binarization is used. The height of the threshold base map obtained by morphological gray filtering in step 2) is increased by a certain iterative value each time as the height threshold for segmenting the DSM to obtain a binary map. The image block is converted into multiple different binary maps, and subsequent steps are then performed to determine whether they are feature points.

[0029]

[0030] F t+1 F(x,y) = t F(x,y)+ΔF

[0031] In the formula, B t (x,y) is the binary image, and the threshold base map is F t (x,y), Z t (x,y) is the original DSM, and ΔF is the iterative value.

[0032] For further optimization, the size of the iterative value needs to be less than the height of the feature, generally set to 0.25m - 0.5m. If the iterative value is too large, it is easy to miss features; if the iterative value is too small, there will be redundant calculations.

[0033] 4) Connected component labeling:

[0034] Perform connected component labeling on the binary image obtained after iteratively binarizing in step 3). Connected component labeling is a basic operation in image processing used to identify and label connected regions formed by adjacent pixels in a binary image. The algorithm starts from the upper left corner of the image and traverses the pixels row by row and column by column. For each pixel traversed, check whether its value is the foreground value 1. If so, perform labeling. At the same time, check its surrounding adjacent pixels. If there are already labeled pixels among the adjacent pixels, the current pixel is labeled as the same connected component. If none of the adjacent pixels are labeled, it means that the current pixel belongs to a new connected component, and a new label is assigned to it.

[0035] For further optimization, this algorithm uses 4-connected labeling, that is, only check the four adjacent pixels above, below, left, and right:

[0036] f 4(x，y) ={f(x + 1, y), f(x - 1, y), f(x, y + 1), f(x, y - 1)}

[0037] f 4(x，y) is the extracted connected component, and f(x, y) is the binary image.

[0038] 5) Threshold filtering:

[0039] For the connected components obtained in step 4), first set a pixel number threshold. Traverse the connected components obtained in the binary image, calculate the number of pixels in each connected component. If the number of pixels is greater than the pixel number threshold, exclude it from subsequent operations. This step is to reduce the computational amount. The set value of the pixel number threshold is generally the number of pixels of the largest feature in this area. Then traverse each pixel in the non-excluded connected components. As long as there are background pixels among the adjacent pixels above, below, left, and right of the pixel, it is identified as an edge pixel, thus achieving the extraction of edge pixels.

[0040] Calculate the average value of the edge pixels of each connected component, compare it with the height difference calculation base map obtained by morphological gray-scale erosion, and set a judgment threshold, preferably 1m - 2m. If the difference between the average value of the edge pixels and the average value of the corresponding pixels in the height difference calculation base map is greater than the judgment threshold, then this connected component is considered as the ground object range, that is, the ground object pixels to be filtered out.

[0041] 6) Surface interpolation:

[0042] Set the ground object pixels to be filtered out obtained in step 5) as null values, then perform interpolation on the area after filtering out the ground object pixels, and then use filtering to smooth the interpolation surface to obtain the ground of the ground object area, that is, the complete DEM.

[0043] 7) Merging process:

[0044] After separately performing steps 2) - 6) on each image block, finally merge the processing results.

[0045] Advantages of the present invention:

[0046] The present invention integrates morphological gray-scale filtering and iterative threshold segmentation techniques, combines a block optimization strategy and a dynamic connected component screening mechanism, and effectively solves problems such as sensitivity to window size, poor terrain adaptability, and a large amount of manual intervention in traditional algorithms. Experiments were verified using the ISPRS Vaihingen dataset. The results show that the average omission rate of ground features is 8.93%, the false detection rate is 3.09%; the average accuracy rate reaches 80.85%, which is 8.07% higher than the average accuracy rate of existing mainstream algorithms; the Kappa coefficient is 75.24%, which is 0.1036 higher than the average of mainstream algorithms in the prior art on average, significantly improving the automation processing efficiency and terrain restoration accuracy under complex terrains, and is applicable to the refined terrain modeling requirements in fields such as geological disaster prediction and water conservancy projects. Brief Description of the Drawings

[0047] The present invention will be further described below in conjunction with the drawings and embodiments.

[0048] Figure 1 It is a flowchart of DEM generation based on the morphological filtering iterative algorithm in Embodiment 1.

[0049] Figure 2 It is a schematic diagram of block processing in Embodiment 1, showing how to perform block processing on large-area images.

[0050] Figure 3 It is a schematic diagram of iterative binarization processing in Embodiment 1, showing how to segment ground features and the ground from the terrain surface by iteratively increasing values.

[0051] Figure 4 It is the position distribution of test areas a, b, and c in Embodiment 1.

[0052] Figure 5 It is a schematic diagram of the manually outlined ranges and processing results of test areas a, b, and c in Embodiment 1, showing the finally generated DEM. Detailed Embodiment

[0053] Embodiment 1:

[0054] A method for generating a digital elevation model (DEM) based on morphological filtering iteration, comprising the following steps:

[0055] Selection of test area and preliminary preparation:

[0056] Collect data. The DSM data is from the ISPRS Vaihingen dataset with a resolution of 25 cm. This area is in a semi-rural and semi-urban state, containing two landforms with sparse and dense buildings, which is quite representative.

[0057] Perform simple denoising on the DSM, mainly by filtering to remove high-frequency and low-frequency noise points to prevent the noise from affecting subsequent algorithm processing of the image.

[0058] 1) Image block division:

[0059] Divide the input image into image blocks according to a certain size, process each image block separately, and finally merge them. Set an overlapping area between blocks. Use the graph cut algorithm for the pixel values in the overlapping area to find the splicing gaps between blocks.

[0060] As Figure 2 shown, in this embodiment, perform block processing on the DSM data to ensure that the size of each block area is larger than the maximum size of the ground objects. The pixel values in the overlapping area are spliced through the graph cut algorithm. Set the block size to 1024×1024 pixels and the overlapping area to 64 pixels wide.

[0061] 2) Morphological gray filtering:

[0062] Perform the opening operation and erosion operation of morphological operations on the image blocks in step 1). Use the structuring element to perform the opening operation on the DSM image to obtain a rough terrain surface, and then use the structuring element to perform the erosion operation to obtain the DSM with the ground object edges eroded. Use the rough terrain surface as the threshold base map for iterative binarization, and use the DSM with the ground object edges eroded as the base map for height difference calculation.

[0063] Use the opening operation and erosion operation to process the divided DSM data. The size of the structuring element for the opening operation is 120×120, and the size of the structuring element for the erosion operation is set to 7×7 to ensure that there are no ground objects in the rough terrain surface and accurate filtering of the DSM ground object edges.

[0064] 3) Iterative binarization:

[0065] Use the method of iterative binarization for the image blocks obtained in step 1). Increase the height value of the threshold base map obtained by the morphological operation in step 2) by a preset iteration value each time as the height threshold for DSM binarization, and convert the image segmentation into multiple different binary images.

[0066] As Figure 3 shown, in this embodiment, increase the height of the threshold base map obtained by morphological gray filtering by a certain iteration value each time as the height threshold for segmenting the DSM to obtain binary images, and convert the image into multiple different binary images. The height value for each iteration is 0.25m to ensure that the ground object range can be successfully detected.

[0067] 4) Connected component labeling:

[0068] Perform connected component labeling on the binary image obtained after iterative binarization in step 3). For each pixel, check whether its value is the foreground value. If it is, perform labeling and also check its surrounding adjacent pixels. If there are labeled pixels among the adjacent pixels, the current pixel is labeled as the same connected component; if none of the adjacent pixels are labeled, it is a new connected component and a new label is assigned to it.

[0069] In this embodiment, connected component labeling is performed on all binary images, and 4-connected labeling is used, that is, only the four adjacent pixels above, below, left, and right are checked to determine whether they belong to the same connected component.

[0070] 5) Threshold filtering:

[0071] For each connected component obtained in step 4), first set a pixel number threshold, and eliminate the connected components with the number of pixels greater than the pixel number threshold; extract the edge pixels of each connected component (if there is a background pixel among the adjacent pixels above, below, left, and right of the pixel, the pixel is recognized as an edge pixel) and calculate the average value. Compare the average value of the edge pixels with the elevation difference calculation base map obtained by morphological gray-scale erosion in step 2). Set a judgment threshold. If the difference between the average value of the edge pixels and the average value of the corresponding pixels in the elevation difference calculation base map is greater than the judgment threshold, the connected component is considered to be the ground object range.

[0072] In this embodiment, the maximum pixel number threshold for ground objects is set to 250000. Traverse the connected components obtained in the binary image, calculate the number of pixels in each connected component. If the number of pixels is greater than the maximum pixel number threshold, it is eliminated from the subsequent calculation, and the connected components with a number less than this are extracted with edge pixels. Calculate the average value of the edge pixels of the connected component, compare it with the average value of the corresponding pixels in the elevation difference calculation base map obtained by morphological gray-scale erosion, and set the threshold to 1m. If the difference is greater than 1m, the connected component is considered to be the ground object range.

[0073] 6) Surface interpolation:

[0074] Set the extracted ground object range to null values, and then use interpolation to fill the null values. The interpolated surface is simply filtered to achieve smoothing, and a digital elevation model (DEM) is obtained.

[0075] 7) Merging process:

[0076] After separately performing steps 2) to 6) on each image block, finally merge the processing results.

[0077] The selected plot and the schematic diagram of the results in this embodiment are as Figure 4 、 5 shown. In Figure 4Among them, the example area is located in the town, with features such as buildings and vegetation of different sizes, and the terrain has a certain undulation; in area a, there are mainly single-family buildings of different sizes connected in a cluster and a small amount of vegetation adjacent to the buildings, with a relatively large slope undulation. In area b, the features are mainly long-strip buildings and dense patches of vegetation, with a very large slope undulation. In area c, the slope is relatively gentle, and the features are mostly single-family buildings and single trees of vegetation, which are relatively scattered. Result verification:

[0078] Compared with the manually marked feature range, taking the accuracy of extracting the feature range as an index, the first type of error, the second type of error, the accuracy rate, and the Kappa coefficient are used to evaluate the algorithm effect, as shown in Table 1.

[0079] Table 1

[0080]

[0081]

[0082] In the test area of the Vaihingen dataset, during the process of generating DEM using this algorithm, compared with the manually marked feature range, the average value of the first type of error (missing detection of features) is 8.93%, the average value of the second type of error (false detection of features) is 3.09%, the average accuracy rate is 80.85%, and the average Kappa coefficient is 75.24%.

[0083] In complex terrain, the method for generating a digital elevation model (DEM) based on morphological filtering iteration shows strong robustness, successfully filtering out most of the features, and the ground curve remains relatively smooth.

[0084] Finally, it should be noted that the above is only used to illustrate the technical solution of the present invention and not to limit it. Although the present invention has been described in detail with reference to the preferred parameter setting scheme, those of ordinary skill in the art should understand that the technical solution of the present invention can be modified or equivalently replaced without departing from the spirit and scope of the technical solution of the present invention.

Claims

1. A method for generating a digital elevation model based on morphological filtering iteration, characterized in that: After collecting and inputting digital surface model (DSM) image data, the following steps are included: 1) Image segmentation: Divide the input DSM image into image blocks; 2) Morphological gray filtering: For the image blocks in step 1), perform an opening operation to obtain a rough terrain surface, then perform an erosion operation to obtain the DSM with the edges of the objects eroded. Use the rough terrain surface as the threshold base map for iterative binarization, and use the DSM with the edges of the objects eroded as the base map for height difference calculation; 3) Iterative binarization: For the image blocks obtained in step 1), use the method of iterative binarization. Increase the height value of the threshold base map obtained from the morphological operation in step 2) by a preset iterative value each time as the height threshold for DSM binarization, and segment the image blocks into multiple different binary images; 4) Connected component labeling: Perform connected component labeling on the binary images obtained after iterative binarization in step 3). For each pixel, check whether its value is the foreground value. If so, perform labeling, and at the same time check its surrounding adjacent pixels; if there are already labeled pixels among the adjacent pixels, the current pixel is labeled as the same connected component as it; if none of the adjacent pixels are labeled, it is a new connected component, and assign a new label to it; 5) Threshold filtering: For each connected component obtained in step 4), first set a pixel number threshold, and exclude the connected components with the number of pixels greater than the pixel number threshold; extract the edge pixels of each remaining connected component and calculate the average value. Compare the average value of the edge pixels with the base map for height difference calculation obtained from the morphological gray erosion in step 2). Set a judgment threshold. If the difference between the average value of the edge pixels and the average value of the corresponding pixels in the base map for height difference calculation is greater than the judgment threshold, consider this connected component as the object range; 6) Surface interpolation: Set the pixels within the object range obtained in step 5) to null values, then perform interpolation on the null value area, and then smooth the interpolation surface using filtering to obtain the digital elevation model (DEM); 7) Merging process: After separately performing steps 2) - 6) on each image block, finally merge the processing results.

2. The method for generating a digital elevation model based on morphological filtering iteration according to claim 1, wherein: In step 1), the size of each image block is set to be larger than the maximum size of the objects in this area; a overlapping area is set between the blocks, and for the pixel values in the overlapping area, use the graph cut algorithm to find the stitching gaps between the blocks.

3. A method for generating a digital elevation model based on morphological filtering iteration according to claim 1, characterized in that: In step 2), in the opening operation, the size of the structural element is set to be larger than the maximum object radius in this area; in the erosion operation, the size of the structural element is set to be no larger than 7×7.

4. A method for generating a digital elevation model based on morphological filtering iteration according to claim 1, characterized in that: In step 3), the size of the preset iterative value is 0.25m - 0.5m.

5. A method for generating a digital elevation model based on morphological filtering iteration according to claim 1, characterized in that: In step 4), the adjacent pixels are the four adjacent pixels above, below, left, and right.

6. A method for generating a digital elevation model based on morphological filtering iteration according to claim 1, characterized in that: In step 5), the method for extracting edge pixels is: as long as there is a background pixel among the adjacent pixels above, below, left, and right of a pixel, then this pixel is identified as an edge pixel.

7. A method for generating a digital elevation model based on morphological filtering iteration according to claim 1, characterized in that: In step 5), the judgment threshold is 1 - 2m.

Images