A method for removing interference objects and filling holes in DEM images

Abstract

The present invention provides a method for removing interference objects and filling holes in DEM images. By performing annotation processing on the DEM images, using the annotated images to perform masking processing on the original images, and for the holes in the masked images, determining the interpolation points by delimiting the boundaries and using the linear interpolation method, further constructing a Delaunay triangulation network through these interpolation points, simultaneously determining the pixel values of the interpolation points, and recursively assigning values to the constructed triangulation network to complete the filling of the image holes. Through this method, for the first time, a hole filling algorithm is introduced in the field of eliminating interference in areas such as vegetation, realizing the correction of the elevation of the ground object area, and finally constructing a high-precision digital elevation model DEM.

Description

Technical Field

[0001] The present invention relates to the field of removing and filling image interference objects, in particular to a method for removing interference objects and filling holes in DEM images. Background Art

[0002] At present, it is relatively mature to produce a digital surface model (DSM) for visible light images using the oblique photogrammetry technology. However, the surface elevation obtained by the oblique photogrammetry technology includes the heights of ground objects. In large-scale engineering designs with a scope of square kilometers or more, the terrain often contains vegetation and other coverings. If contour lines are directly generated using the digital surface model DSM, elevation anomaly areas will appear and the true terrain conditions cannot be reflected. Therefore, in order to restore the ground elevation, it is usually necessary to filter out the heights of ground objects such as vegetation. And in the actual process of filtering out ground objects such as vegetation, it is easy to have phenomena such as incomplete filtering or poor filtering effect, which causes difficulties for subsequent filling work. In addition, in the field of filling holes in actual images, especially in the filtering of ground objects such as vegetation, a complete and feasible filling scheme has not been proposed for actual operation. Summary of the Invention

[0003] The present invention provides a method for removing interference objects and filling holes in DEM images, and for the first time proposes a hole filling algorithm in the field of eliminating ground objects such as vegetation, so as to realize the correction of the elevation of the ground object area.

[0004] The technical means adopted by the present invention are as follows:

[0005] A method for removing interference objects and filling holes in DEM images, which removes interference objects in the DEM image and fills holes through the following steps:

[0006] Step 1: Perform annotation processing on the DEM image for interference objects, and perform masking processing according to the annotated image to obtain a DEM image with holes;

[0007] Step 2: Extract the minimum rectangular boundary of the area where the holes are located in the DEM image obtained in Step 1, and appropriately expand it to obtain an index file;

[0008] Step 3: Taking the expanded rectangular area in Step 2 as the boundary, determine the positions of interpolation points according to the known points within the rectangular area by the linear interpolation method, and construct a Delaunay triangular mesh surface that satisfies the Delaunay criterion through the determined interpolation points;

[0009] Step 4: While constructing the Delaunay triangular mesh surface in Step 3, determine the pixel values of the interpolation points according to the pixel values of the known points within the rectangular area by the linear interpolation method;

[0010] Step 5: Using the positions and pixel values of the interpolation points obtained in Steps 3 and 4, assign values to the entire Delaunay triangular mesh surface using a recursive algorithm to obtain the filled DEM image.

[0011] Preferably, the outer expansion distance of the minimum rectangular boundary in Step 2 is 2 pixels.

[0012] Preferably, after obtaining the expanded rectangular area in Step 2, further record the positions of the upper left corner and the lower right corner of this rectangular area, and then obtain the index file in Step 2.

[0013] Preferably, the order of determining the interpolation points and constructing the triangular mesh surface in Steps 3 and 4 is: from left to right and from top to bottom along the expanded rectangular area.

[0014] Preferably, the marking process for the interfering objects in Step 1 includes: adding a class name and a class label value to the interfering objects in the DEM image.

[0015] Preferably, the masking process in Step 1 includes masking the masked image according to the masking image of the marked interfering objects, and obtaining the DEM image with holes after removing the interfering objects through the masking process analysis tool.

[0016] Preferably, in the masking process,

[0017] The masking image is: a vector masking map with interfering object markings obtained by vectorizing the DEM image in Step 1 and performing interfering object marking processing on the vectorized image;

[0018] The masked image is: a digital surface model matching the DEM in Step 1.

[0019] The present invention has the following advantages compared with the existing technologies:

[0020] This method uses the combination of the label classification map and the digital elevation model for masking processing in ground object coverage processing for the first time. Compared with traditional ground object coverage processing, the masking effect of this method is better; and this method proposes a multi-hole DEM image filling algorithm based on irregular triangular meshes for the first time in the production of large-scale topographic maps (1:500 - 1:2000), realizes the correction of elevation anomalies caused by trees, etc., and finally constructs a high-precision digital elevation model DEM. Description of the Drawings

[0021] Figure 1 It is a schematic flow chart of the ground elevation extraction method based on visible light images and deep learning algorithms.

[0022] Figure 2It is a top view of a surveyed area including a vegetation area in an embodiment.

[0023] Figure 3 is Figure 2 the corresponding contour map.

[0024] Figure 4 It is the corresponding surveyed area range and distribution map of phase control points for obtaining photos in an embodiment.

[0025] Figure 5 is for Figure 4 the orthophoto DOM obtained by processing the photos in the surveyed area.

[0026] Figure 6 is for Figure 4 the digital surface model DSM obtained by processing the photos in the surveyed area.

[0027] Figure 7 It is the selected sample area dom13N after format conversion and normalization processing of the orthophoto DOM.

[0028] Figure 8 It is a Figure 7 sample area slope13N corresponding to the position of the sample area selected after format conversion and normalization processing of the digital surface model DSM.

[0029] Figure 9 It is the vector image obtained after annotation processing of the multi-channel image.

[0030] Figure 10 is Figure 9 the corresponding raster data image.

[0031] Figure 11 It is the sample data set obtained after shearing and enhancement processing of the multi-channel image and the class label map.

[0032] Figure 12 The data (sample) set to be predicted obtained after shearing and enhancement processing of the four-channel image to be predicted.

[0033] Figure 13 It is the prediction result vector map corresponding to the image to be predicted.

[0034] Figure 14 It is the DEM image containing holes obtained after masking processing according to the prediction result vector map.

[0035] Figure 15 It is the superimposed map of the contour line corresponding to the DEM image obtained after hole filling and the contour line generated from the original digital surface model DSM to be measured. Specific implementation mode

[0036] Specifically, in combination with the attached drawings of the specificationFigures 1-15 , the following specific solutions are given:

[0037] Combined with a topographic mapping project of a certain ash field, this method is described. The survey area of the project is about 12 square kilometers, and the mapping scale is 1:500. In this project, the Pegasus D2000 drone is used in the fieldwork, equipped with a SONY a6000 camera, and integrated with a POS and IMU system with RTK function. The sensor size of the SONY a6000 camera is 23.5X15.5mm, with an effective pixel of 24 million and a focal length of 25mm.

[0038] In this project, the survey area is crisscrossed with ravines, such as Figure 2 shown, which is a top view of the survey area of this embodiment including the vegetation area. As shown in the framed area, there are many tall trees scattered in the gully. If the digital surface model DSM is directly used to generate the contour lines as Figure 3 shown, it can be seen that there are multiple elevation anomaly areas (black masses), which cannot reflect the true terrain conditions and bring errors to the cut and fill calculation.

[0039] Therefore, for the above problems, according to the actual operation steps, as Figure 1 shown, this method proposes a specific embodiment:

[0040] First step, use the drone equipped with an optical camera to obtain a large number of centimeter-level spatial resolution photos in a selected area. The selected area mainly uses vegetation-sparse areas at the square kilometer level. It is necessary to synchronously obtain the internal and external orientation element information of the photos and synchronously arrange photo control points that meet the specification requirements.

[0041] Specifically, the aerial photography time of this embodiment project is May 17, 2020, and the weather is clear. The flight altitude is about 300 meters, the photo overlap degree in the flight direction ≥ 70%, preferably 80%, and the photo overlap degree in the flight side direction ≥ 60%, preferably 65%. A total of 3619 photos are taken in this project, the average spatial resolution of the photos is 0.05m, the range of a single photo is about 280m * 212m, and a total of 29 photo control points are arranged in the survey area, as Figure 4 shown, which is the distribution map of the survey area and photo control points of this project. Among them, Figure 4 in (a) is the distribution map of the survey area and photo control points, Figure 4 in (b) is the photo control point.

[0042] Second step, use the drone data processing software to input the photos, photo control points and internal orientation element information obtained in the first step into the above drone data processing software to obtain an orthophoto map (DOM) and a digital surface model (DSM) that match the survey area. The orthophoto map DOM and digital surface model DSM obtained in this step are the basic images for subsequent processing. Among them, the drone data processing software uses Context Capture software. AsFigure 5 and Figure 6 are respectively the orthophoto DOM and digital surface model DSM (in tif format) corresponding to the survey area processed by using this software.

[0043] In the third step, convert the digital surface model DSM obtained in the second step into the corresponding slope map. Generally, the obtained orthophoto DOM is an 8-bit unsigned integer of RGB three channels, and the digital surface model DSM is a 32-bit floating point type. And since both the orthophoto DOM and the digital surface model DSM are relatively sensitive to vegetation, for example, vegetation in the orthophoto DOM is manifested as spectral and texture features, and the digital surface model DSM is manifested as abrupt changes in elevation gradient and texture features. Here, in this method, the slope is calculated pixel by pixel for the digital surface model DSM to obtain the slope map corresponding to this selected area. The slope calculation considers the 8 neighborhoods of the central pixel. The following Table 1 shows the schematic diagram of the positions of the central pixel e and the surrounding 8 pixels:

[0044]

[0045]

[0046] Table 1

[0047] Assume that a, b, c, d, e, f, g, h, i are respectively elevation values, and x_cellsize and y_cellsize are the common unit distances in the x direction and the y direction. To calculate the slope of the central pixel e, it is necessary to calculate the change rate of pixel e in the x direction (Formula (1)) and the change rate of pixel e in the y direction (Formula (2)) respectively:

[0048] [dz / dx] = ((c + 2f + i) - (a + 2d + g)) / (8 * x_cellsize) (1)

[0049] [dz / dy] = ((g + 2h + i) - (a + 2b + c)) / (8 * y_cellsize) (2)

[0050] Furthermore, the slope value at pixel e is obtained (Formula (3)):

[0051]

[0052] Moreover, by performing format conversion and normalization on the two types of data, a 4-channel floating-point data image map containing RGB channels and slope channels is obtained. This method is the first to simultaneously use two types of data, the orthophoto DOM and the slope map, in the deep learning algorithm for vegetation classification. Combining vegetation classification files such as class label maps, deep learning extraction for vegetation classification is carried out. Of course, this method can also be applied to the classification of other ground objects or interference objects. In this embodiment, the vegetation classification in mapping is mainly described, and other channel images can also be introduced on this basis to integrate into a multi-channel image convenient for subsequent vegetation classification and extraction.

[0053] Among them, format conversion means converting both the orthophoto and the slope map into the same floating-point data, and integrating the three RGB channels and the slope channel of the orthophoto to obtain a 4-channel map; the normalization process is to uniformly take values between 0 and 1. Specifically, a slope map slope.tif (value range 0.0000 - 90.0000) is generated according to dsm.tif, and it is normalized to obtain a single-channel floating-point data slopeN.tif (value range 0.0000 - 1.0000). The orthophoto DOM is in the format of 8-bit unsigned integer for three RGB channels (value range 0 - 255), and each channel data is normalized to obtain a 3-channel floating-point data domN.tif (value range 0.0000 - 1.0000). A region is selected from each of the 3-channel floating-point data domN.tif corresponding to the orthophoto DOM and the single-channel floating-point data slopeN.tif corresponding to the slope map as the sample area, which are respectively the sample area dom13N.tif as shown in Figure 7 and the sample area slope13N.tif as shown in Figure 8 . And the data is further spliced to obtain a 4-channel data sample.tif (value range 0.0000 - 1.0000) containing R, G, B, and slope, keeping the projection and the image rows and columns the same as the original Figure 1 .

[0054] Step 4: Perform label annotation processing on ground objects such as vegetation for the four-channel image obtained in Step 3 (or use the orthophoto obtained in Step 2). Specifically, use remote sensing data processing software to interpret the image. The remote sensing data processing software selects ArcGIS software, and carefully annotates the corresponding vegetation areas to generate vector polygons matching the image, add fields of class names and class label values to its attribute values (such as background 0, vegetation 1, and more classes can be extended according to the type of ground objects), and save it as a vector file (label.shp) as shown in Figure 9 , and further use ArcGIS software to convert the vector file label.shp into the same size (same resolution and number of rows and columns) and projection method as the four-channel image asFigure 10 The raster file shown, namely the class label map (label.tif).

[0055] Step 5: Perform image clipping on the class label map label.tif obtained in Step 4 and the four-channel image sample.tif generated in Step 3, that is, clip the classification label map obtained in Step 4 and the multi-channel image obtained in Step 3 into sample maps of the same size, and at the same time perform sample enhancement processing (such as rotation, blurring, adding noise, etc.). One feature of the enhancement processing is to enhance the display of the parts to be recognized that are marked, which is convenient for better learning and recognition; the enhancement processing obtains a sample dataset of four-channel maps and corresponding class label maps of the same size and meeting the size of an integer multiple of 32, such as Figure 11 The sample dataset obtained after clipping and enhancing the four-channel image and the class label map, with a total of 6,000 images, and the image size is 512 * 512 pixels each.

[0056] Step 6: Randomly divide the sample dataset generated in Step 5 into two parts, one part for model training, that is, the training dataset, and one part for model verification, that is, the verification dataset.

[0057] In this specific embodiment, it is preferably to divide the sample dataset generated in Step 5 into two parts according to a ratio of 9:1. The part with a ratio of 9 is used for model training, and the part with a ratio of 1 is used for model verification. And save the file name indexes of the training dataset and the verification dataset to a file.

[0058] Step 7: Perform multi-channel recognition and ground object classification recognition training on the training dataset in Step 6 through the UNET model. UNET can obtain efficient classification results with a small number of training samples. It adjusts parameters according to task requirements, including the number of data channels, the number of categories, the number of iterations, the learning rate, etc. In this embodiment, the structure of the backbone network adopts the VGG16 model, uses the Relu function as the activation function, and uses the cross-entropy and softmax functions to calculate the loss. Use the training dataset generated in Step 6 and the initial model file (transferable learning) as inputs, and perform 400 rounds of training to finally implement a UNET vegetation classification model that can perform 4-channel data training and prediction of R, G, B, and slope, and generate a series of initial model files (*.pth).

[0059] In the eighth step, use the validation dataset in the sixth step and the initial model file generated in the seventh step as inputs for validation training. Calculate the mean intersection over union (MIOU) of the predicted values and the true values for each channel and each ground object category recognition. This is one of the reference factors for model update. Modify the parameters of the initial model based on this MIOU. Preferably, the calculated MIOU is excellent. Otherwise, modify the model parameters, iterate and optimize the model, and update the model file (*.pth) to obtain the final prediction recognition model (M). The calculation steps of the MIOU are as follows:

[0060] (1) Calculate the confusion matrix;

[0061] (2) Calculate the intersection over union IOU_i for each channel and ground object category. The formula is:

[0062] IOU_i = TP_i / (TP_i + FN_i + FP_i)

[0063] Where:

[0064] TP_i is the true positive value, that is, the intersection of the true value and the predicted value;

[0065] FN_i is the false negative value, that is, the part of the true value after removing the intersection of the true value and the predicted value;

[0066] FP_i is the false positive value, that is, the part of the predicted value after removing the intersection of the true value and the predicted value;

[0067] (3) Average the IOU of each class to obtain the mean intersection over union MIOU.

[0068] In the ninth step, perform the processing of the second and third steps on the image to be predicted to obtain a four-channel image to be predicted corresponding to the image to be predicted, which includes the RGB three channels and the slope channel.

[0069] In the tenth step, then perform the shearing and enhancement processing of the fifth step on the four-channel image to be predicted obtained in the ninth step to obtain the image to be predicted data (sample) set corresponding to the image to be predicted as shown. The size of each sample is the same as the size of the sample obtained after the shearing and enhancement processing of the four-channel image and the class label map obtained in the fifth step. Also, the obtained image to be predicted data file needs to retain the projection and geographic coordinate information of the original image to be measured to facilitate positioning and stitching. The naming rule of the data file is "filename_original total number of rows_original total number of columns_starting row number_starting column number.tif". Figure 12 shown, and the size of each sample is the same as the size of the sample obtained after the shearing and enhancement processing of the four-channel image and the class label map obtained in the fifth step. And, the obtained image to be predicted data file needs to retain the projection and geographic coordinate information of the original image to be measured to facilitate positioning and stitching. The naming rule of the data file is "filename_original total number of rows_original total number of columns_starting row number_starting column number.tif".

[0070] In the eleventh step, send the image to be predicted data set generated in the tenth step and the recognition model file obtained in the eighth step into the final prediction recognition model for prediction to obtain the vegetation prediction result data set.

[0071] Step 12: Mosaic (i.e., splice) the prediction results obtained in Step 11 to obtain prediction results of the same size and projection as the original four-channel image to be predicted in Step 9.

[0072] Step 13: Convert the raster of the prediction results obtained in Step 12 into a vector, such as Figure 13 to obtain a vector map of the prediction results corresponding to the image to be predicted, from which the boundary range of the vegetation in the circle can be obtained.

[0073] Step 14: Overlay the vector file of the prediction results obtained in Step 13 on the four-channel image to be predicted in Step 9 for interactive editing to add, delete, and modify missed and misdetected areas, further improving the vegetation recognition rate, and obtaining a ground object classification result map with full coverage for multi-channel information recognition and ground object classification recognition such as vegetation.

[0074] Step 15: Use the repaired classification results obtained in Step 14 as a mask to perform masking on the digital surface model DSM processed in the UAV data processing software corresponding to the image to be predicted, to obtain a Figure 14 DEM image containing several holes (vegetation areas) as shown.

[0075] Step 16: Based on the DEM image with holes obtained in Step 15, extract the minimum rectangular boundary of each vegetation (or ground object) area in the DEM image, and expand it appropriately (by 2 pixels), record its upper left and lower right positions, and obtain the corresponding index file.

[0076] Step 17: According to the index file obtained in Step 16, fill the holes one by one from left to right and from top to bottom. For each sub-region with holes, use the linear interpolation method for three-dimensional spatial interpolation. That is, use the linear interpolation method to determine the interpolation points based on the known points such as the rectangular boundary determined in Step 16, and further, based on triangles, first find 3 points around the interpolation point to form a triangle according to the Delaunay method. The interpolation point is inside the triangle, and thus a Delaunay triangular mesh surface that satisfies the Delaunay criterion is constructed through the determined interpolation points. Use this method to obtain a continuous and smooth surface, that is, the elevation surface. And, during the process of constructing the elevation surface, the algorithm retains the original pixel values and fills in the missing pixel values, that is, determines the pixel values of the interpolation points according to the pixel values of the known points in the rectangular region by the linear interpolation method, and recursively assigns values to the entire Delaunay triangular mesh surface according to the position and pixel values of the interpolation points, to obtain a Figure 15The smooth and continuous DEM data image of the research area to be measured is shown. In this figure, the interrupted lines are the contour lines obtained by using this method, and the solid lines are the contour lines directly generated using the digital surface model DSM. This method first proposed a multi-hole DEM image filling algorithm, filling the gap in this field.

[0077] This method is applicable to the extraction of DEM in areas with sparse vegetation, and can especially be used for the generation of large-scale topographic maps (1:500 - 1:2000) in areas of square kilometers. Through this method, the true ground elevation value can be kept unchanged, vegetation can be filtered out, and the ground elevation of the corresponding area can be filled. It has a high degree of automation and can greatly improve the work efficiency of industry insiders.

Claims

1. A method for removing interference objects and filling holes in DEM images, characterized in that, Remove the interference objects in the DEM image and fill the holes through the following steps: Step 1: Mark the interference objects in the DEM image, and perform masking processing based on the marked image to obtain a DEM image with holes; Step 2: Extract the minimum rectangular boundary of the area where the holes are located in the DEM image obtained in Step 1, and appropriately expand it to obtain an index file; Step 3: Taking the expanded rectangular area in Step 2 as the boundary, determine the positions of the interpolation points by the known points within the rectangular area according to the linear interpolation method, and construct a Delaunay triangular mesh surface that meets the Delaunay criterion through the determined interpolation points; Step 4: While constructing the Delaunay triangular mesh surface in Step 3, determine the pixel values of the interpolation points by the pixel values of the known points within the rectangular area according to the linear interpolation method; Step 5: Using the positions and pixel values of the interpolation points obtained in Steps 3 and 4, assign values to the entire Delaunay triangular mesh surface using a recursive algorithm to obtain the filled DEM image; Among them, the marking process of the interference objects in the DEM image in Step 1 includes: adding class names and class label values to the interference objects in the DEM image; The masking process in Step 1 includes masking the masked image according to the masked image of the marked interference objects, and processing through the masking processing analysis tool to obtain a DEM image with holes after removing the interference objects; In the masking process, the masked image is: a vector masked map with interference object markings obtained by vectorizing the DEM image in Step 1 and performing interference object marking processing on the vectorized image; the masked image is: a digital surface model matching the DEM in Step 1.

2. The method for removing interference objects and filling holes in a DEM image according to claim 1, characterized in that, The expansion distance of the minimum rectangular boundary in Step 2 is 2 pixels.

3. A method for removing interference objects and filling holes in a DEM image according to claim 2, characterized in that, After obtaining the expanded rectangular area in Step 2, further record the positions of the upper left corner and the lower right corner of the rectangular area, and then obtain the index file in Step 2.

4. A method for removing interference objects and filling holes in a DEM image according to claim 1, characterized in that, The order of determining the interpolation points and constructing the triangular mesh surface in Steps 3 and 4 is: from left to right and from top to bottom along the expanded rectangular area.

Images