A method for generating DSM by optimizing grid using satellite stereo image pairs

Abstract

The present invention relates to a method for generating a DSM by optimizing a grid using a satellite stereo image pair. First, a rough original three-dimensional grid is generated from a multi-scene linear array satellite stereo image pair and used as input. Then, based on the RPC model representing the relationship between object points and image points, reprojection errors and Laplacian gradients are calculated, etc. Two energy terms, namely the photometric consistency energy term and the gradient smoothing energy term, are respectively constructed with these two calculated quantities as the main indicators to complete the construction of the energy function. The three-dimensional grid is continuously subdivided with an adaptive threshold. A method of discrete first and then optimization is used for the optimization solution process. This method avoids the problems of the size adaptability of the dense matching window and the implicit problem of the disparity consistency within the window in the traditional generation method, and has higher accuracy.

Description

Technical Field

[0001] The present invention belongs to the field of 3D modeling, and particularly relates to a method for generating a DSM by optimizing a grid using a satellite stereo image pair. Background Art

[0002] In the production process of map products, the Digital Surface Model (DSM) is an important pre - data source for tasks such as the extraction of the Digital Elevation Model (DEM), the generation of contour lines, the extraction and reconstruction of buildings, the production of true orthophotos, and the update of geographic information. The existing methods for generating DSM mainly include Light Detection and Ranging (LiDAR) and photogrammetry.

[0003] The method of directly obtaining DSM by LiDAR has a high acquisition efficiency and is not restricted by sunlight and weather. However, the directly obtained 3D point cloud has "bottlenecks" in post - processing technologies such as the removal of interference points, the collection of break lines, and the interpolation of the building surface for maintaining discontinuity.

[0004] The key technology for obtaining DSM by photogrammetry depends on the dense matching of images. Window - correlation - like methods are generally used, but they face problems such as the adaptability of the relevant window size and the implicit restriction of the parallax consistency within the window. Although the semi - global matching (SGM) method has achieved an effect comparable to LiDAR, its computational cost and memory overhead are large, and it depends on the maximum disparity search range and the number of pixels of the image. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a method for generating a DSM by optimizing a grid using a satellite stereo image pair. Through this method, it is possible to drive the generation of DSM by means of grid optimization when only a rough grid is input. The present invention provides a DSM generation method that is different from the traditional dense matching method of photogrammetry and has better effects.

[0006] To achieve the above object, based on the grid optimization method of linear array satellite images, the problem of DSM generation is transformed into a 3D grid optimization problem. Based on the RPC model of the linear array satellite stereo image pair, the problem of optimizing the driving of 3D grid vertices based on general images is transformed into the problem of directly driving the optimization of 3D grid vertices using linear array satellite images, and the optimized 3D grid (DSM) is obtained by constructing an energy function for optimization and solution.

[0007] The technical solution of the present invention is as follows:

[0008] First, a rough original three-dimensional grid is generated from a linear array satellite stereo image pair. Based on the RPC model representing the relationship between object points and image points, the reprojection error and Laplacian gradient are calculated. Using these two calculated quantities as the main indicators, the photometric consistency energy term and the gradient smoothing energy term of the energy function are respectively constructed to complete the construction of the energy function. The three-dimensional grid is subdivided with an adaptive threshold, and the optimization solution process is carried out using the method of discrete first and then optimization. At the same time, the grid is continuously subdivided to continuously optimize the original grid and drive the vertices to move to positions close to those corresponding to the actual research area. After the iterative solution ends, its x, y, and z coordinates are close to the three-dimensional coordinates of the corresponding points in the real three-dimensional world, and the final optimization result is the DSM.

[0009] This method is different from traditional photogrammetric dense matching. It directly uses the linear array satellite image pair to drive the vertices of the three-dimensional grid to achieve grid optimization, and the optimization result is the generated DSM.

[0010] Specifically, the method of the present invention includes the following steps:

[0011] (1) First, obtain the linear array satellite stereo image pair of the target area and its RPC (Rational Polynomial Coefficients) parameters, and obtain a rough three-dimensional grid (with longitude, latitude, and elevation as coordinates) through the grid construction method;

[0012] (2) Based on the initial three-dimensional grid and the linear array satellite stereo image pair, calculate the photometric consistency energy term and the gradient smoothing energy term of the energy function, and construct the energy function;

[0013] E(S) = λ photo E photo (S) + λ smooth E smooth (S)

[0014] S is the surface of the three-dimensional grid;

[0015] E(S) is the energy function of the three-dimensional grid;

[0016] E photo (S) is the photometric consistency energy term;

[0017] E smooth (S) is the energy term for gradient smoothing of the three-dimensional grid;

[0018] λ photo is the weight of the photometric consistency energy term; λ smooth is the weight of the gradient smoothing energy term;

[0019] λ photo and λ smoothThe weight value defaults to 0.5 and can be changed according to the final result of iterative optimization. When this value becomes smaller, the corresponding λ smooth is large, and the gradient smoothing term controls the smoothness of the established 3D mesh. Therefore, it will make the 3D mesh more tend to be smooth. The weight value can be controlled to adjust the smoothness of the final mesh. The photometric consistency energy term is calculated with ZNCC (Zero-normalized cross-correlation) as the main index, and the gradient smoothing energy term is calculated with the first-order and second-order Laplacian gradients as the main indexes;

[0020] (3) Subdivide the entire 3D mesh with an adaptive threshold to make the mesh details richer and closer to the real ground surface;

[0021] Use the area of the triangular mesh of the 3D mesh projected onto the 2D image as the judgment value (actually the sum of pixels included in the projection area). The threshold is default set to 24. When the projection value is greater than this threshold, the triangular mesh is subdivided. The subdivision judgment is performed for each iteration. When the iteration number is reached, the subdivision stops simultaneously.

[0022] (4) Solve by the gradient descent method, update iteratively, drive the synchronous movement of the triangular mesh vertices, and obtain the optimal solution, that is, optimize the triangular mesh to generate DSM.

[0023] The present invention proposes to directly optimize the rough original grid using linear array satellite images, convert the corresponding relationship between object points and image points through the RPC model, and then calculate ZNCC to characterize the photometric consistency energy term in the energy function. This term is used to measure the offset degree of a point when it is re-projected from one image to another through the 3D mesh; at the same time, calculate the gradient of the 3D mesh surface, which is used to measure the smoothness of the 3D mesh surface; then apply regularization constraints to balance the first two terms (photometric term and smooth term); finally, optimize and solve the constructed energy function (adjust the positions of the triangular mesh vertices constituting the three-dimensional mesh by the gradient descent method).

[0024] The method of the present invention directly optimizes the rough grid constructed by linear array satellite images without using the commonly used dense matching method in photogrammetry. Directly optimize the initial rough grid, and a 3D mesh (DSM) close to the real situation of the ground surface can be optimized from an almost blank grid; it provides a new method and idea for DSM production, avoiding the problems of the adaptability of the window size of dense matching and the implicit problem of the disparity consistency within the window in the traditional generation method, and has higher accuracy. Description of the Drawings

[0025] Figure 1 is the rough original grid input. 1 is the contour of a certain building, and 2 is the hole existing in the original grid.

[0026] Figure 2 Schematic diagram of driving the movement of triangular mesh vertices by mesh optimization. 1 is the vertex of the triangle that constitutes the grid, 2 is the basic unit that constitutes the mesh, and 3 is the movement of the vertex.

[0027] Figure 3 Schematic diagram of reprojection error.

[0028] Figure 4 Schematic diagram of three-dimensional mesh subdivision. 1, 2, and 3 are the midpoints of the three sides, and 4, 5, and 6 are the four triangles generated after the triangular mesh is subdivided.

[0029] Figure 5 Schematic diagram after three-dimensional mesh subdivision. Compared with Figure 2 the original mesh, the number of triangular meshes is more.

[0030] Figure 6 is Figure 1 the DSM map generated after mesh optimization.

[0031] Figure 7 Flow chart for generating DSM of satellite stereo image pairs based on mesh optimization. Detailed implementation mode

[0032] The present invention relates to a method for generating DSM by using satellite stereo image pairs for mesh optimization. Based on the RPC model of linear array satellite stereo image pairs, the problem of optimizing the driving of three-dimensional grid vertices based on general images is converted into the problem of directly using linear array satellite stereo image pairs to optimize the driving of three-dimensional grid vertices, and the optimized three-dimensional grid (DSM) is obtained by constructing an energy function for optimization and solution. First, a rough original three-dimensional grid is generated according to the linear array satellite stereo image pairs; then, the corresponding relationship between object points and image points is determined based on the RPC model, and the reprojection error and gradient of the image are calculated to represent the components of the energy function; the three-dimensional grid is continuously subdivided with an adaptive threshold; and an optimization and solution process is carried out by using the method of discrete first and then optimization.

[0033] Based on the above principle, the implementation mode of the present invention will be described in detail with reference to the accompanying drawings. This embodiment is a method for generating DSM of satellite stereo image pairs based on mesh optimization.

[0034] 1. Obtain the initial mesh

[0035] 1.1 Obtain multiple scenes of linear array satellite images of a certain area, and construct a network through general photogrammetry methods to obtain a very rough three-dimensional grid, which can be rough enough for the purpose of serving as the initial plane for mesh optimization. When constructing the initial three-dimensional grid, it can be generated by common dense matching methods (such as SGM) in photogrammetry or open source projects such as ColMap, VisualSFM, and CloudCompare. Only a sufficiently rough three-dimensional grid is required as a processing platform for carrying out.

[0036] As shown in this embodiment Figure 1 is an initial plane (sufficiently rough) generated by using the SGM algorithm in the photogrammetric dense matching algorithm with a high - resolution satellite stereo image pair. Among them, the total number of vertices of this initial grid is 4443, and the number of faces forming the grid is 8671. The entire grid is very rough, with the building outline faintly visible and obvious hole conditions existing.

[0037] 1.2 Pre - processing of the initial grid. Pre - process the initial rough three - dimensional grid, such as projection transformation, etc. Pre - process according to the coordinate system of the constructed initial grid. For example, when the input initial grid is in the longitude - latitude - elevation coordinate system, it is necessary to perform projection transformation to convert it into projected coordinates. Offset the coordinates before processing and restore the coordinates after processing (add or subtract according to the RPC parameters during offset). The purpose is to make the three - dimensional coordinate units consistent for easy display and ensure the accuracy after input.

[0038] In this embodiment, the input is the initial three - dimensional grid generated by a high - resolution satellite stereo image pair. First, perform x and y coordinate offsets, and the offset parameters are: - 81.65295044719999850713, 30.32576204769999961286 (determined according to the latoffset and longoffset parameters in the RPC file). The unit of the initial grid is longitude - latitude and elevation. Through projection transformation, convert the longitude - latitude of UTM 17N into the projected coordinates of UTM 17N, so that the units of x, y, and elevation z all become meters.

[0039] 2. Construct the energy function

[0040] 2.1 Construct the photometric - consistency energy term. Calculate the reprojection error with ZNCC as the main index. The purpose is to measure the deviation degree when the points of one scene of the image are reprojected to another scene of the image through the three - dimensional grid. For example Figure 3 , when the vertices of the reference image are reprojected to another image (Reprojected image) through the grid surface S, there is a reprojection error. Thus, construct the photometric - consistency energy term. Its formula is as follows:

[0041]

[0042] E photo (s) is the photometric - consistency energy term;

[0043] h(I, J)(x) is a decreasing function (opposite to the normalized cross - correlation ZNCC) that measures the photometric consistency between image I and image J;

[0044] represents image I jReproject through the mesh S to I i superior.

[0045] Indicates the valid area of ​​reprojection;

[0046] According to the above formula, the photometric consistency energy term is constructed, mainly expressed in ZNCC, and is calculated as follows:

[0047]

[0048] Among them, f(x, y) is the original image, t(x, y) is the template image, n is the number of pixels (elements) in the template, σ is the standard deviation of the sample, and μ is the mean of the sample.

[0049] In this embodiment, ZNCC is calculated by calculating the reprojection error of the right image reprojected to the left image through the three-dimensional grid, and h(I, J)(x) is represented by ZNCC to constitute the photometric consistency energy term of the energy function.

[0050] 2.2 Construct the gradient smoothing energy term. The specific implementation is to calculate the Laplace gradient, the purpose is to ensure that the three-dimensional grid is relatively smooth and penalize the curvature. The construction formula is as follows:

[0051]

[0052] k1, k2 are the principal curvatures of the surface where the point is located.

[0053] 2.3 Reasonable allocation of λ photo and λ smooth The energy function is now constructed.

[0054] In this embodiment, λ photo and λ smooth The weights of are set to 0.5, 0.5 respectively, and their sum is 1. photo and λ smooth Respectively control the realism and smoothness of the three-dimensional grid. smooth When it increases, the final 3D mesh edges are smoother. Its value is adjusted according to the final result, and the weight value is continuously adjusted to achieve a balance between realism and smoothness.

[0055] 3. Adaptive mesh subdivision

[0056] Subdivide the 3D mesh according to the judgment conditions. When the number of pixels corresponding to the image reprojection onto another image exceeds the threshold and the texture complexity exceeds the threshold (the sum of the number of pixels contained in the projection area and the texture complexity of the projection area), divide the triangle. The determination of the subdivision threshold is determined by the projection relationship between each triangular face constituting the 3D grid and the image. Calculate the sum of the number of pixels and the texture complexity within the projection area of each triangular face on the image, and set it as the threshold. If it is greater than this threshold, subdivision is performed. The division principle is as Figure 4 shown. Take the midpoint of each edge as a new vertex, and divide the original triangle into four new triangles. The purpose is to achieve mesh subdivision and make the details more perfect.

[0057] In this embodiment, as Figure 5 shown, the parameter of the subdivision threshold is set to 24, and this threshold can be set manually. With the change of the parameter setting, the fineness of the final mesh is also different. The larger the threshold, the coarser the mesh details; on the contrary, the richer the mesh details.

[0058] 4. Optimization of the energy function

[0059] 4.1 Optimization solution of the energy function. The initial value of the iteration is the initial 3D mesh. In actual use, the gradient descent method is used for iteration to obtain the optimal solution. Each time the iteration calculates the offset of each vertex of the triangle, and continuously iterates until the vertices of the triangular mesh move to the optimal position close to the true situation of the ground surface. As Figure 5 shown, the entire 3D mesh is composed of thousands of triangular meshes. When using the gradient descent method for solution, the initial gradient is the set of all vertices in the mesh, and each vertex is composed of the coordinates in the x, y, and z directions. During the solution process, the values in the x, y, and z directions change iteratively, causing the vertex coordinates to move continuously. Stop when the set number of iterations is reached, and the result at the stop is the optimal result. The schematic diagram of the vertex movement is as Figure 2 .

[0060] In this embodiment, the iteration parameter is set to 255.05, and the step size is set to 1.5. The optimized result after movement is as Figure 6 . Compared with the original input mesh Figure 1 , the level of detail is significantly richer, the building outlines are clearly visible, and the holes in the original mesh are also filled.

[0061] 4.2 Limit the step size and move synchronously. During the iterative solution process, to prevent the movement of a certain vertex in the triangular mesh from being too large, it is necessary to limit the movement of each vertex. Take the average side length of the three-dimensional mesh as the threshold to constrain the movement of the driving vertex, so that each vertex moves synchronously. Calculate the side length of each triangular mesh that makes up the entire three-dimensional mesh as the basis for judging the limit step size. When the movement in the x, y, and z directions is greater than the average side length of each triangular mesh, limit the movement (gradient change) in the x, y, and z directions to half of the average side length of the three-dimensional mesh. This enables the three vertices of the triangular mesh to move synchronously, avoiding the situation where the solution result of the gradient descent method is poor due to the excessive movement of a certain vertex.

[0062] In this embodiment, when the movement of a certain vertex is too large, under the condition of a certain number of iterations, it may be impossible to optimize to the optimal position. For example, Figure 5 as shown, the bulge in the circled area. Considering the reason, it may be that the movement of the vertices of several triangular meshes is too large compared to other vertices, so the "bulge" and "depression" situations occur. Therefore, it is necessary to limit the movement of the vertices.

Claims

1. A method for generating a DSM by optimizing a grid using a satellite stereo image pair, characterized in that, First, generate a rough original three-dimensional grid based on the linear array satellite stereo image pair. Represent the relationship between object points and image points according to the RPC model, calculate the reprojection error and the Laplacian gradient, and use these two calculation quantities as the main indicators to respectively construct the photometric consistency energy term and the gradient smoothing energy term of the energy function, thus completing the construction of the energy function; Subdivide the three-dimensional grid with an adaptive threshold, and adopt the method of discretization first and then optimization to carry out the optimization solution process; At the same time, continuously subdivide the grid to continuously optimize the original grid, and drive the vertices to move to the position closest to the actual situation. The final optimization result is the DSM.

2. The method according to claim 1, characterized in that, The method includes the following steps: (1) First, obtain the linear array satellite stereo image pair of the target area and its RPC parameters, and obtain a rough three-dimensional grid through the grid construction method; (2) Based on the initial three-dimensional grid and the linear array satellite stereo image pair, calculate the photometric consistency energy term and the gradient smoothing energy term of the energy function, and construct the energy function; ； S is the surface of the three-dimensional grid; is the energy function of the three-dimensional grid; is the energy term for photometric consistency; The energy term for three-dimensional grid gradient smoothing; is the weight of the photometric consistency energy term; is the weight of the gradient smoothing energy term; Calculate the photometric consistency energy term with ZNCC as the main indicator, and calculate the gradient smoothing energy term with the first and second order Laplacian gradients as the main indicators; (3) Subdivide the entire three-dimensional grid with an adaptive threshold to make the grid details richer and closer to the real ground surface; (4) Solve by the gradient descent method, iterate and update, drive the triangular mesh vertices to move synchronously, and obtain the optimal solution, that is, the optimized triangular mesh, and generate the DSM.

3. The method according to claim 2, wherein: ； is a decreasing function that measures the photometric consistency between and the photographic film ; Presentation photograph Through the grid Reprojected onto Above; Indicates the valid region for reprojection; Construct the photometric consistency energy term according to the above formula, mainly represented by ZNCC, and calculate as follows: ； Among them, is the original image, is the template image, and n is the number of pixels in the template; By calculating the reprojection error of the right image reprojected through a three-dimensional grid onto the left image, the ZNCC is calculated, and the ZNCC is used to characterize , which constitutes the photometric consistency energy term of the energy function; ； is the principal curvature of the surface where the point is located; and The weight value is defaulted to 0.5.

Images