A method for generating digital elevation models based on morphological filtering iteration

Abstract

The present invention discloses a method for generating a digital elevation model based on iterative morphological filtering, comprising the following steps: 1) image segmentation; 2) morphological grayscale filtering; 3) iterative binarization; 4) connected domain labeling; 5) threshold filtering; 6) surface interpolation; and 7) merging. The present invention integrates morphological grayscale filtering with iterative threshold segmentation techniques, combines a segmentation optimization strategy with a dynamic connected domain screening mechanism, and effectively addresses the issues of traditional algorithms such as window size sensitivity, poor terrain adaptability, and the need for multiple parameter settings. This method significantly improves the automated processing efficiency and terrain restoration accuracy of generating a digital elevation model (DEM) from a digital surface model (DSM) in complex terrain.

Description

Technical Field

[0001] The present invention relates to a digital elevation model (DEM) generation technology, belonging to the field of geographic information processing technology. Specifically, it provides a digital elevation model generation method based on morphological filtering iteration, which can achieve efficient and accurate digital surface model object filtering and DEM generation. Background Art

[0002] Digital elevation models (DEMs) are important data models that reflect the surface topography and are widely used in fields such as geological disaster prediction, water conservancy project construction, and environmental protection. Currently, DEM generation typically relies on digital surface model (DSM) data. However, DSMs contain objects such as buildings and vegetation, which need to be removed to obtain an accurate DEM.

[0003] In the field of digital terrain model processing, the traditional method for converting digital surface models (DSM) into digital elevation models (DEM) is manual editing, which uses empirical judgment to screen and remove ground objects. This method consumes a lot of labor and time, and requires experienced workers to ensure accuracy.

[0004] With the advancement of automatic computer image processing, many DSM processing algorithms have been applied, and a slope-based filtering algorithm has emerged. This algorithm assumes that the slope value between ground points and non-ground points is greater than the slope value between ground points, thereby filtering out ground features in the DSM. Its advantages are high speed and efficiency. However, the single fixed slope threshold makes it perform poorly in mountainous areas with complex terrain.

[0005] There is also morphological filtering, which applies morphological opening operations to compare the height changes of corresponding points before and after the opening operation. Points exceeding the threshold are considered to be ground objects. However, the algorithm depends on the size of the sliding window. If it is too small, it cannot filter out large areas of ground objects, and if it is too large, it will filter out undulating terrain.

[0006] Later, to address the shortcomings of the morphological algorithm, progressive morphological filtering was proposed. Sliding windows of different sizes were used to perform morphological operations. Finally, the results of different sizes were weighted to obtain the final filtering result. This better solved the problem of difficulty in determining the window size, but the top of the area with a larger slope was still over-filtered.

[0007] The progressive irregular triangulated network filtering algorithm selects initial ground seed points to construct a triangulated network, and determines whether unclassified points are added to the triangulated network through conditions. It iterates continuously until no new points are added to the triangulated network, and uses the points that constitute the triangulated network to generate a DEM to achieve the effect of removing ground objects. This algorithm has high filtering accuracy, but it depends on the selection of initial seed points and requires low iteration efficiency.

[0008] In recent years, cloth simulation filtering has been proposed based on the principles of physical simulation. This algorithm flips all points upside down. 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 to cope with different terrain conditions, and it performs poorly in areas with complex terrain types. Summary of the Invention

[0009] To address the shortcomings of the existing algorithms mentioned above in the DEM generation process, the present invention provides an iterative algorithm based on morphological filtering that can effectively remove objects from the DSM and generate a high-precision digital elevation model (DEM). This method uses a morphological algorithm to obtain a terrain surface, iteratively increases the terrain surface height value to cut the surface into a binary image sequence. The height difference between the edge extracted from the binary image sequence and the pixel of the original image is used to determine the range of the object. The object is then filtered out and interpolated, thereby achieving the filtering of the objects in the DSM and the reconstruction of the DEM surface.

[0010] The object of the present invention is achieved like this:

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

[0012] 1) Image segmentation:

[0013] The input digital surface model (DSM) image is divided into image blocks of a certain size (each image block is processed separately in the subsequent steps, and the processing results are finally merged). As an optimization, overlapping areas are set between blocks, and the pixel values ​​in the overlapping areas are used using 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 even hundreds of millions of pixels, loading and processing them block by block can not only improve memory efficiency and avoid the problem of insufficient device memory, but also reduce computational complexity and ease the difficulty of overall processing.

[0015] When setting the block size for segmentation, the invention sets each image block to a size larger than the maximum size of the terrain features in the area. Specifically, the number of vertical and horizontal pixels in the image block exceeds the maximum number of vertical and horizontal pixels of the feature outline. Furthermore, to ensure smooth transitions between image blocks, a certain overlap area is set between the blocks. A graph cut algorithm is used to calculate the pixel values ​​in the overlapping area to find the appropriate gap between the blocks. This graph cut algorithm optimizes the inter-block splicing and ensures terrain continuity.

[0016] Graph Cut is an image segmentation algorithm based on the minimum energy path. It treats the image as a graph and uses the minimum cut to find the optimal splicing boundary, thereby achieving 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 area, 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 cutting seam to minimize ∑ for all s,t on t he cut M(s,t,A,B), the pixels on both sides of the cutting seam are taken from the A image and the B image respectively.

[0019] 2) Morphological grayscale filtering:

[0020] The image block in step 1) is subjected to morphological opening and erosion operations, and the DSM image is opened using a structuring element to obtain a rough terrain surface. The DSM with the edges of the objects eroded is then eroded using the structuring element. The rough terrain surface is used as the threshold base map for iterative binarization, and the DSM with the edges of the objects eroded is used as the base map for height difference calculation.

[0021] Morphological grayscale filtering is an extension of mathematical morphology and is mainly used to process grayscale images. The present invention uses the opening and erosion operations of morphological grayscale 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, the DSM image is opened using a suitable structuring element to obtain a rough terrain surface. During the opening operation, the structuring element size must be larger than the maximum feature radius in the area. Otherwise, features will appear raised within the rough terrain, affecting subsequent iterative binarization. Therefore, a larger structuring element should be used whenever possible. An erosion operation is then performed to obtain a DSM with the feature edges eroded away. In this step, the structuring element should not be set too large; set the structuring element (filter) size to 7×7 or smaller. The two results serve as the threshold basemap for iterative binarization and the basemap for height difference calculation, respectively.

[0026] 3) Iterative binarization:

[0027] To filter out the ground objects, we need to find a binary height threshold that can separate the ground objects from the ground. The image is binarized using this threshold to separate the ground points from the non-ground points.

[0028] The present invention uses an iterative binarization method to apply the iterative binarization method to the image block obtained in step 1), and increases the height of the threshold base map obtained by morphological grayscale filtering in step 2) by a certain iteration value each time, which is used as the height threshold of the binary image obtained by segmenting the DSM, and converts the image block into multiple different binary images, and then performs subsequent steps to determine whether it is a ground feature point.

[0029]

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

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

[0032] For further optimization, the iteration value needs to be smaller than the height of the feature, generally set to 0.25m-0.5m. If the iteration value is too large, features may be missed, while if the iteration value is too small, redundant calculations may occur.

[0033] 4) Connected domain labeling:

[0034] Connected domain labeling is performed on the binary image obtained after iterative binarization in step 3). Connected domain 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 traversed pixel, its value is checked to see if it is a foreground value of 1. If so, it is labeled. At the same time, the surrounding adjacent pixels are checked. If there are labeled adjacent pixels, the current pixel is labeled as having the same connected domain as it. If none of the adjacent pixels are labeled, it means that the current pixel belongs to a new connected domain and is assigned a new label.

[0035] For further optimization, this algorithm uses 4-connected labeling, that is, only checking 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 connected component extracted, and f(x, y) is a binary image.

[0038] 5)Threshold filtering:

[0039] For the connected domain obtained in step 4), a pixel count threshold is first set. The connected domains obtained in the binary image are traversed, and the number of pixels in each connected domain is calculated. If the number of pixels is greater than the pixel count threshold, it is removed from subsequent calculations. This step is to reduce the amount of calculation. The pixel count threshold is generally set to the number of pixels of the largest feature in the area. Then, each pixel in the connected domain that has not been removed is traversed. If there are background pixels in the adjacent pixels above, below, left, or right of the pixel, it is identified as an edge pixel, thus achieving edge pixel extraction.

[0040] Calculate the average value of the edge pixels of each connected domain and compare it with the height difference calculation base map obtained by morphological grayscale erosion. Set a judgment threshold, preferably 1m to 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, the connected domain is considered to be the object range, that is, the object pixel to be filtered out.

[0041] 6) Surface interpolation:

[0042] The ground object pixels to be filtered out obtained in step 5) are set to null values, and then the area where the ground object pixels are filtered out is interpolated. After that, the interpolation surface is smoothed by filtering to obtain the ground surface of the ground object area, that is, the complete DEM.

[0043] 7) Merge processing:

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

[0045] Beneficial effects of the present invention:

[0046] This method integrates morphological grayscale filtering with iterative threshold segmentation techniques, combined with a block optimization strategy and a dynamic connected domain screening mechanism, effectively addressing issues such as window size sensitivity, poor terrain adaptability, and frequent manual intervention in traditional algorithms. Experimental validation using the ISPRS Vaihingen dataset revealed an average missed detection rate of 8.93% and a false detection rate of 3.09%. The average accuracy reached 80.85%, an 8.07% improvement over the average accuracy of existing mainstream algorithms. The Kappa coefficient reached 75.24%, an average improvement of 0.1036 over the average accuracy of mainstream algorithms in existing technologies. This method significantly improves automated processing efficiency and terrain restoration accuracy in complex terrain, making it suitable for the refined terrain modeling needs of fields such as geological disaster prediction and water conservancy projects. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The present invention will be further described below with reference to the accompanying drawings and examples.

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

[0049] Figure 2 This is a schematic diagram of block processing in Example 1, showing how to process a large area image in blocks.

[0050] Figure 3 This is a schematic diagram of the iterative binarization process in Example 1, showing how to segment the terrain surface from the ground by iteratively increasing the value.

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

[0052] Figure 5 This is a schematic diagram of the manually outlined range and processing results of test areas a, b, and c in Example 1, showing the final generated DEM. DETAILED DESCRIPTION

[0053] Example 1:

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

[0055] Selection of test area and preliminary preparation:

[0056] The DSM data were collected from the ISPRS Vaihingen dataset with a resolution of 25 cm. The area is semi-rural and semi-urban, with both sparsely built and densely built landforms, which is quite representative.

[0057] The DSM is subjected to simple denoising, mainly to remove high-frequency and low-frequency noise through filtering to prevent the subsequent algorithm processing of the noisy image.

[0058] 1) Image segmentation:

[0059] The input image is divided into image blocks of a certain size, each image block is processed separately, and finally merged. The overlapping area is set between the blocks, and the pixel values ​​of the overlapping area are used using the graph cut algorithm to find the splicing gaps between the blocks;

[0060] like Figure 2 As shown, in this embodiment, the DSM data is divided into blocks to ensure that the size of each block is larger than the maximum size of the ground feature, and the pixel values ​​of the overlapping area are spliced ​​using the graph cut algorithm. The block size is set to 1024×1024 pixels, and the overlapping area is 64 pixels wide.

[0061] 2) Morphological grayscale filtering:

[0062] The image block in step 1) is subjected to morphological opening and erosion operations, and the DSM image is opened using a structuring element to obtain a rough terrain surface. The DSM with the edges of the objects eroded is then eroded using the structuring element. The rough terrain surface is used as the threshold base map for iterative binarization, and the DSM with the edges of the objects eroded is used as the base map for height difference calculation.

[0063] The block-based DSM data were processed using opening and erosion operations. The size of the structural element for the opening operation was 120 × 120, and the size of the structural element for the erosion operation was set to 7 × 7 to ensure that there were no objects in the rough terrain surface and that the edges of the DSM objects were accurately filtered out.

[0064] 3) Iterative binarization:

[0065] The image blocks obtained in step 1) are subjected to an iterative binarization method. The threshold base map height value obtained by the morphological operation in step 2) is increased by a preset iteration value each time as the height threshold of the DSM binarization, and the image segmentation is converted into multiple different binary images.

[0066] like Figure 3 As shown, in this embodiment, the height of the threshold base map obtained by morphological grayscale filtering is increased by a certain iteration value each time, which is used as the height threshold of the binary image obtained by segmenting the DSM, and the image is converted into multiple different binary images. The height value of each iteration is 0.25m, ensuring that the range of the ground object can be successfully detected.

[0067] 4) Connected domain labeling:

[0068] The binary image obtained after iterative binarization in step 3) is labeled with connected domains. For each pixel, check whether its value is a foreground value. If so, mark it and check its surrounding adjacent pixels at the same time. If there are marked pixels among the adjacent pixels, the current pixel is marked as the same connected domain as it. If no adjacent pixels are marked, it is a new connected domain and a new label is assigned to it.

[0069] In this embodiment, connected domains are marked for all binary images using 4-connected marking, that is, only the four adjacent pixels (upper, lower, left, and right) are checked to determine whether they belong to the same connected domain.

[0070] 5)Threshold filtering:

[0071] For each connected domain obtained in step 4), first set a pixel number threshold and eliminate the connected domains with pixel numbers greater than the pixel number threshold; extract the edge pixels of each connected domain (as long as there are background pixels in the adjacent pixels above, below, left, and right of the pixel, the pixel is identified as an edge pixel) and calculate the average value. The average value of the edge pixels is compared with the height difference calculation base map obtained by morphological grayscale erosion in step 2), and a judgment threshold is set. 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, the connected domain is considered to be the range of the ground object.

[0072] In this example, the maximum pixel count threshold for ground objects is set to 250,000. The connected domains obtained in the binary map are traversed, and the number of pixels in each connected domain is calculated. If the number of pixels is greater than the maximum pixel count threshold, it is removed from subsequent calculations. Connected domains with a number less than this threshold are extracted as edge pixels. The average value of the edge pixels in the connected domain is calculated and compared with the average value of the corresponding pixels in the height difference calculation base map obtained by morphological grayscale erosion. The threshold is set to 1 meter. If the difference is greater than 1 meter, the connected domain is considered to be the ground object range.

[0073] 6) Surface interpolation:

[0074] The extracted ground feature range is set to null value, and then the null value is filled by interpolation. The interpolated ground surface is smoothed by simple filtering to obtain the digital elevation model (DEM).

[0075] 7) Merge processing:

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

[0077] The plots selected in this embodiment and the result diagram are as follows: Figure 4 、 5 As shown. Figure 4In the example, the area is located in a town, with buildings of varying sizes and vegetation, and the terrain is somewhat undulating. Area a is mainly distributed with single buildings of varying sizes and a small amount of vegetation adjacent to the buildings, with a large slope. Area b is mainly composed of long strips of buildings and dense vegetation, with a large slope. Area c has a relatively gentle slope, and the features are mostly single buildings and single plants, with a relatively scattered distribution. Results verification:

[0078] Compared with the manually marked feature range, the accuracy of the extracted feature range is used as an indicator, and the first type error, second type error, accuracy rate, and 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, when using this algorithm to generate DEM, compared with the range of manually marked objects, the average first-type error (missed objects) is 8.93%, the average second-type error (falsely detected objects) is 3.09%, the average accuracy is 80.85%, and the average Kappa coefficient is 75.24%.

[0083] In complex terrain, the digital elevation model (DEM) generation method based on morphological filtering iteration shows strong robustness, successfully filters out most of the ground objects, and keeps the ground curve relatively smooth.

[0084] Finally, it should be noted that the above is only used to illustrate the technical solution of the present invention and is not limiting. Although the present invention is described in detail with reference to the preferred parameter setting solution, ordinary technicians in this field should understand that the technical solution of the present invention can be modified or replaced by equivalents 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 by: After collecting and inputting the digital surface model (DSM) image data, the following steps are included: 1) Image segmentation: Divide the input DSM image into image blocks; 2) Morphological grayscale filtering: For the image block in step 1), an opening operation is performed to obtain a rough terrain surface, and then an erosion operation is performed to obtain a DSM with the edges of the objects eroded away. The rough terrain surface is used as the threshold base map for iterative binarization, and the DSM with the edges of the objects eroded away is used as the base map for height difference calculation; 3) Iterative binarization: The image block obtained in step 1) is subjected to an iterative binarization method. The threshold base map height value obtained by the morphological operation in step 2) is increased by a preset iteration value each time, which is used as the height threshold of the DSM binarization. The image block is segmented and converted into multiple different binary images. 4) Connected domain labeling: The binary image obtained after iterative binarization in step 3) is labeled with connected domains. For each pixel, check whether its value is a foreground value. If so, label it. At the same time, check its neighboring pixels. If there are marked pixels among the neighboring pixels, the current pixel is marked as having the same connected domain as it. If no neighboring pixels are marked, it is a new connected domain and is assigned a new label. 5) Threshold filtering: For each connected domain obtained in step 4), first set a pixel number threshold and eliminate connected domains with pixel numbers greater than the pixel number threshold; extract the edge pixels of each remaining connected domain and calculate the average value; compare the average value of the edge pixels with the height difference calculation base map obtained by morphological grayscale erosion in step 2), and 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 height difference calculation base map is greater than the judgment threshold, the connected domain is considered to be a ground feature range; 6) Surface interpolation: Set the pixels within the feature range obtained in step 5) to null values, then interpolate the null value areas, and then use filtering to smooth the interpolation surface to obtain the digital elevation model (DEM); 7) Merging process: After performing steps 2) to 6) on each image block separately, the processing results are finally merged.

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 ground object in the area; overlapping areas are set between blocks, and the pixel values ​​in the overlapping areas are used using a graph cut algorithm to find the splicing gaps between blocks.

3. The method for generating a digital elevation model based on morphological filtering iteration according to claim 1, wherein: Step 2) In the opening operation, the size of the structure element is set to be larger than the radius of the maximum feature in the area; in the erosion operation, the size of the structure element is set to be no larger than 7×7.

4. The method for generating a digital elevation model based on morphological filtering iteration according to claim 1, wherein: Step 3) The preset iteration value is 0.25m-0.5m.

5. The method for generating a digital elevation model based on morphological filtering iteration according to claim 1, wherein: In step 4), the adjacent pixels are four adjacent pixels at the top, bottom, left, and right.

6. The method for generating a digital elevation model based on morphological filtering iteration according to claim 1, wherein: Step 5) The edge pixel extraction method is: as long as there are background pixels in the adjacent pixels above, below, left, and right of the pixel, the pixel is identified as an edge pixel.

7. The method for generating a digital elevation model based on morphological filtering iteration according to claim 1, wherein: The judgment threshold in step 5) is 1 to 2 m.

Images