Research on Building DSM Fusion Method Based on Adaptive Spline and Target Characteristic Guidance

Abstract

In order to adapt to the actual scene of a stereo satellite observing the same area sequentially and improve the accuracy of the target-oriented 3D reconstruction, this paper proposed a building DSM fusion update method based on adaptive splines and target characteristic guidance. This method analyzed the target characteristics of surface building targets to explore their intrinsic geometric structure information, established a nonlinear fusion method guided by the target characteristics to achieve the effective fusion of multiple DSMs on the basis of maintaining the target structural characteristics, and supported the online updating of DSM to ensure the needs of practical engineering applications. This paper presented a DSM fusion method for surface building targets and finally conducted DSM fusion experiments using typical urban area images of different scenes. The experimental results showed that the proposed method can effectively constrain and improve the DSM of buildings, and the integrity of the overall construction of the target 3D model structure was significantly improved, indicating that this paper provides an effective and efficient DSM constraint method for buildings.

1. Introduction

Three-dimensional (3D) reconstruction is a research hotspot in the application of computer-aided engineering technology for product and engineering design [1,2,3]. Three-dimensional (3D) reconstruction is a new technology that integrates knowledge in the fields of machine vision, image processing, and computer cartography and has characteristics that cannot be compared with 2D pictures. Traditional remote sensing 3D reconstruction methods usually aim at large-area and large-scale scene reconstruction and tend to realize the model construction of the overall structure of the target from a global perspective, commonly without focusing on the local structure of the target, and it is difficult to ensure the structural integrity and accuracy of the target with complex structure [4,5,6].

In the field of target 3D reconstruction, the main data sources at this stage include digital surface model (DSM) [7], digital elevation model (DEM), digital terrain model (DTM), and airborne LiDAR data [8] generated from aerial imagery and satellite high-resolution image interpretation. In the past decades, researchers have proposed many methods to implement target 3D model construction based on DSM or LiDAR data, which can be summarized into three categories: model-driven, data-driven, and hybrid-driven methods [9,10,11].

Model-driven methods start with the construction of a library of predefined primitives containing the basic model of the target, based on which the target 3D model construction is achieved by matching the structure between the predefined model and the input data [12]. Model-driven methods are commonly applied to the 3D reconstruction process of models whose targets are projected as rectangles in the 2D plane. Vosselman et al. [13] used geometric invariant moment features for primitive model matching and estimated the parameters of the model primitives to achieve 3D model construction. Subsequently, Brenner et al. [14] proposed a splitting rule for target shape and structure and realized the construction of a three-dimensional model of the target with a more complex structure. The model construction method proposed by Henn et al. [15] used the SVM algorithm to statistically classify the point clouds of buildings within the rectangularly divided area of the target bottom surface and the RANSAC algorithm to fit the model parameters to achieve the estimation of the target model structure parameters to obtain the final target construction results. Lafarge et al. [16] combined DSM data and the extracted 2D contours of the target for 3D model reconstruction, which is based on the concept of model splitting and stitching, using the Gibbs model to control the combination mode between each primitive model of the target, using Markov Chain Monte Carlo (MCMC) to find the optimal combination mode between model primitives, and estimating the structural parameters of each primitive model, so as to realize the 3D reconstruction of the overall structural model of the target. Huang et al. [17] proposed a generalized model construction method that treats the construction of the target as a top-down process, predefines the parametric description and combination method of the primitive models in the model library, and uses the Reversible Jump Markov Chain Monte Carlo (RJMCMC) technique combined with the structural transformation kernel function to achieve 3D model construction of targets with different structures. However, it is difficult to construct targets with complex structures based on model-driven target 3D reconstruction methods. Due to the complex and variable target geometric structure, it is difficult to parametrically describe the 3D model of targets with different degrees of complexity by a unified method, which will affect the completeness and accuracy of the overall structural model construction of the target without the constraints of the global structure.

In recent years, data-driven model construction methods have been widely used in the field of 3D reconstruction due to their higher flexibility [18,19,20]. The data-driven model construction method first uses segmentation methods to obtain different adjacent subsurfaces of the top surface of the target and improves the regularization of the overall structure of the target by adjusting the structure of adjacent points and lines between the subsurfaces to obtain the final model construction results. Verma et al. [21] proposed a data-driven model construction method, which introduced a seed region selection method based on the similarity of local features based on the traditional segmentation method of region growth to achieve the segmentation of the top surface of the target building by effectively controlling the range of region growth through the constraints of the structural features of local regions. Finally, the results of the 3D reconstruction of the model are obtained by combining the boundary points of the target top surface and the connection topology of the ridge line features. Chauve et al. [22] proposed a region growing method based on topological constraints on the target top surface for segmentation, which added structural constraints to the topological connections of the target top surface to the traditional region growing algorithm, thus enabling more accurate acquisition of information on each subplane included in the top surface of the target, which improved the completeness of the target model construction. Sampath et al. [23] used the fuzzy K-mean clustering technique to segment the target top surface and obtain the boundary and ridge features of each subplane and adjusted the connection between each vertex, as well as the line features according to the predefined normalized structure criterion to achieve the model construction of the overall structure of the target. Kim et al. [24] used the level-sets approach for model construction, which firstly constructed an energy function based on multitemporal level sets and obtained the segmentation results of the target top surface according to the minimization process of solving the energy function and then obtained the contours and ridges of the target top surface according to the output of the zero-level energy function and finally obtained the final model 3D reconstruction results according to the boundaries of the adjacent subsurfaces of the top surface, the ridges, and the topological connection relationships between the vertices. The limitation of the data-driven approach is that the quality of model construction is too dependent on the segmentation results of the target top surface. When there are certain noise or singularities in the data source, the segmentation for the target top surface usually has the phenomenon of over-segmentation or undersegmentation, resulting in the final construction results losing local structural features, and thus affecting the structural integrity of the target construction model.

In order to further improve the stereo model construction, a hybrid stereo model construction method combining data-driven and model-driven has been widely used [25,26], which mainly combines the advantages of model-driven and data-driven methods and firstly constructs a library of well-formed building primitive models based on the model-driven concept and predefines the parameterization of different types of model primitive description methods. Secondly, the data-driven method is used to detect structural features such as boundaries and ridges of the target, and finally, the 3D reconstruction of the target is realized by the constraints of the connection rules of these basic structural elements. Xiong et al. [27] proposed a dictionary library based on the structured cartographic description of the model, which stores the structured cartographic drawings for each type of structure, and the dictionary library has representative error subdrawings and corresponding structure correction subdrawings, based on which the structural integrity and accuracy of the model 3D reconstruction are further improved by the techniques of cartographic matching and structure correction.

In summary, traditional 3D reconstruction methods are aimed at building large regional terrain or urban scenes, and commonly neglect the fine construction of the local structure of the target. In the process of 3D reconstruction for the target, there is a higher demand for the completeness and accuracy of the overall construction of the target 3D model structure, and more focus on the high-precision construction of the local structure of the target, which are the challenges faced by the traditional 3D reconstruction methods.

The remainder of this paper is organized as follows. Section 2 details the specific workflow of research, including building geometry characteristics, adaptive spline-based DSM fusion guided by building geometric characteristics. Section 3 is a comparative analysis of experimental results. Discussions and conclusions are presented in Section 4 and Section 5, respectively.

2. Methodology

In this paper, we first analyzed the target characteristics of the surface building target to explore the geometric structure information inherent in the building, established a nonlinear fusion method guided by the target characteristics to achieve the effective fusion of multiple DSMs on the basis of maintaining the target structural characteristics, proposed the DSM fusion method for the surface building target, and finally conducted DSM fusion experiments using typical urban area images of different scenes.

2.1. Building Geometric Structure Characteristics

A typical building is composed of beam–slab–column elements, and its roof is mostly composed of colorful steel sheets or modified bitumen. In this paper, we considered exactly these types of buildings. From the geometric point of view, the roof shape of a typical building is a polyhedral structure, a representative example of which is shown in Figure 1 (derived from reference [17], titled “A generative statistical approach to automatic 3D building roof reconstruction from laser scanning data”).

According to Figure 1 (derived from reference [17], titled “A generative statistical approach to automatic 3D building roof reconstruction from laser scanning data”), it can be found that typical building roof shapes are diverse and contain flat roofs, pitched roofs, double-slope roofs, four-slope roofs, and multiwave folding slab roofs. Although the specific manifestations of roof shapes are diverse, they can be considered to consist of roof panels with different inclination angles. For this reason, it is necessary to construct a unified mathematical description of building roof characteristics to accommodate the different roof geometric structures of typical buildings.

For building roofs with different morphologies, this paper constructed adaptive 2D spline functions to describe building geometry universally. For the 2D spline function, each roof panel of a building can be considered as a linear spline with specific parameters in 3D space, the flat roof of a building can be considered as a 0th-order 2D spline, the pitched roof with different inclinations can be represented as a 1st-order spline with different parameters, and the roof of the whole building can be considered as a complex and diverse roof shape with linear splines connected by ridge lines with different parameters, as shown in Figure 2.

As can be seen from Figure 2, each roof panel of a four-slope building can be considered as an inclined surface of different angles in 3D space, and the top flat roof panel can be represented as a 0th-order spline, while the four sloping roof panels on the side can be represented as 1st-order splines with different parameters. Connecting the spline functions through the ridge line forms the roof geometry of the building, which means that the spline functions can accurately describe a typical building roof.

For typical flat roofs, pitched roofs and typical buildings with flat or pitched roof panels such as multiwave folding roofs all can be accurately described using the spline function, and the conclusion can be summarized as follows.

Theorem: For a rectangular interval defined in the subspace $\left(x,y\right)$ of the right-angle coordinate system $\left(x,y,h\right)$, there exists a set of linear spline functions $S_i:I_i\to ℝ1\le i\le N$ that can accurately fit the roof geometry of typical flat roofs, pitched roofs, and multiwave folding roofs with flat or pitched roof panels under a specific interval division ${I_i|}_{i=1}^N$, and satisfying $I_i\cap I_j=\varnothing ,i\ne j$, ${\displaystyle \underset{i}{\cup }{I}_i}=I$ conditions.

Note that the spline function is a segmental function. Only by dividing ${I_i|}_{i=1}^N$ in the appropriate interval can the specific building roof structure be accurately described. The spline function of each region constitutes each roof panel of the building, the intersection line between the intervals corresponds to the ridge line of the building, and for buildings with different roofs, the location and direction of the ridge line are not fixed. In order to characterize the different directions and positions of the roof ridge lines of different structural buildings, the building area needs to be gridded at different scales to approximate the roof ridge lines of different buildings. The gridding of the building area at different scales is shown in Figure 3.

It can be seen from Figure 3 that the different scales of meshing of the rectangular area of the building can be adapted to buildings with different roof structures. Specifically, a flat-roofed building can be represented as a 0th-order spline defined by a rectangular grid, a pitched-roofed building can be represented as two 1st-order splines defined by two rectangular grids (separated by the projection of the ridge line at $\left(x,y\right)$), and a four-slope roof can be represented as five 1st-order splines defined by four trapezoids with one rectangular grid (separated by the projection of the ridge line at $\left(x,y\right)$). Therefore, the spline function only needs a suitable grid division so that the grid line and the ridge line in the $\left(x,y\right)$ projection line coincide to ensure that the spline function can accurately describe the corresponding structure of the building. The gridding of scales I, II, and III in the figure can be adapted to flat-roofed buildings, pitched-roofed buildings, and four-slope-roofed buildings. It can be inferred that the denser the gridding, the more accurate the spline function can be in describing the complex roof shape of the building.

2.2. Analysis of Fitting Error of Adaptive Spline to Building Target

In practice, building targets have a wide variety of roof shapes, which means that it is difficult to prespecify the exact grid to achieve an accurate description of arbitrarily shaped roofs. In order to complete the meshing of the building area more effectively, the fitting error of the spline function to the typical building roof shape is analyzed.

For any given building roof structure, that is, any n^th-order differentiable function $F:I\to ℝ$ defined in the rectangular interval $I$ in the subspace $(x,y)$ of the right-angle coordinate system $(x,y,h)$, for any given positive number $\epsilon$, there exists a sufficiently dense interval division ${I_i|}_{i=1}^N$ that satisfies $I_i\cap I_j=\varnothing ,i\ne j$, ${\displaystyle \underset{i}{\cup }{I}_i}=I$ so that the linear spline function $S_i:I_i\to ℝ1\le i\le N$ interpolates any noncollinear three points of the building roof panel, and the elevation error at any location of the spline function and the building roof structure is less than $\epsilon$.

Proof: First consider the spline function $S_1:I_1\to ℝ$ in the interval $I_1$. If it is able to interpolate any noncollinear three points $\left(x_1,y_1,h_1\right)$, $\left(x_2,y_2,h_2\right)$ and $\left(x_3,y_3,h_3\right)$ of the building roof, the parameters of this linear spline $a_0+a_{1}x+b_{1}y=h,\left(x,y\right)\in I_1$ can be calculated by the following equation.

$$\left[\begin{array}{l}1,x_1,y_1\\ 1,x_2,y_2\\ 1,x_3,y_3\end{array}\right]\times \left[\begin{array}{l}{a}_0\\ a_1\\ b_1\end{array}\right]=\left[\begin{array}{l}{h}_1\\ h_2\\ h_3\end{array}\right]$$

(1)

Since a unique plane (i.e., the plane described by the linear splines) can be determined through any noncollinear three points, the location of any point on the plane is an affine combination of the three points, which can be expressed as the following equation.

$$S_1\left({\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i,y_i\right]}\right)={\displaystyle \sum _i{l}_i\cdot F\left(x_i,y_i\right)}{\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i,y_i\right]}\in I_1$$

(2)

To calculate the elevation difference between the spline function and the building roof function, Equation (3) is expanded using Taylor’s formula.

$$\begin{array}{l}{S}_1\left({\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i,y_i\right]}\right)={\displaystyle \sum _i{l}_i\cdot (F\left(x_0,y_0\right)+\nabla F{\left(x_0,y_0\right)}^T\left[x_i-x_0;y_i-y_0\right]+}\\ 0.5\times {\left[x_i-x_0;y_i-y_0\right]}^{T}G\left(F\left(x_c,y_c\right)\right)\left[x_i-x_0;y_i-y_0\right])\\ F\left({\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i,y_i\right]}\right)=F\left(x_0,y_0\right)+\nabla F{\left(x_0,y_0\right)}^T\left({\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i-x_0;y_i-y_0\right]}\right)+\\ 0.5\times {\left({\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i-x_0;y_i-y_0\right]}\right)}^{T}G\left(F\left(x_c,y_c\right)\right)\left({\displaystyle \sum _{i=1}^3{l}_i\cdot \left[x_i-x_0;y_i-y_0\right]}\right)\end{array}$$

(3)

Notice that ${\left[x_i-x_0;y_i-y_0\right]|}_{i=1}^3\le l_{len}$, $l_{len}$ is determined by the size of the spline function grid, so the upper limit of the deviation between the actual elevation of the building roof at any location and the elevation fitted by the spline function is $l_{len}{}^2\times C$, C is defined in the following equation and can be approximated as a constant in general. Therefore, the elevation deviation is inversely proportional to the square of $l_{len}$, that is, the denser the mesh, the smaller the maximum elevation deviation. For any given elevation error $\epsilon$, when the grid is sufficiently dense, that is, $l_{len}\le \sqrt{\frac{\epsilon }{C}}$, the maximum elevation error is less than $\epsilon$. This means that a sufficiently dense grid division enables the spline function to fit an arbitrary function accurately.

$$\begin{array}{l}C=\frac{1}{l_{len}{}^2}{\left[\begin{array}{l}{x}_1-x_0,y_1-y_0\\ x_2-x_0,y_2-y_0\\ x_3-x_0,y_3-y_0\end{array}\right]}^T(\left[i_1;i_2;i_3\right]G\left(F\left(x_c,y_c\right)\right){\left[i_1;i_2;i_3\right]}^T+\\ diag\left(i_1,i_2,i_3\right)G\left(F\left(x_c,y_c\right)\right)diag{\left(i_1,i_2,i_3\right)}^T)\left[\begin{array}{l}{x}_1-x_0,y_1-y_0\\ x_2-x_0,y_2-y_0\\ x_3-x_0,y_3-y_0\end{array}\right]\end{array}$$

(4)

2.3. DSM Fusion Guided by Building Geometric Characteristics Based on Adaptive Spline

The spline function with sufficiently dense meshing can accurately describe buildings with arbitrary roof shapes, including not only typical flat, pitched, and multiwave roofs, but also unconventional buildings such as curved and spherical roofs. From an architectural point of view, the design of the roof structure of a building has to consider the waterproof performance and material shaping, which makes the roof of a building mostly flat and pitched. When the roofs are flat and pitched, the spline function with dense meshing can effectively fit the roofs, but the local deformation of the building generated by the inaccurate elevation in DSM will also be fitted by the spline function. In order to maintain the building roof structure as much as possible and avoid the introduction of inaccurate elevation information in the DSM, this paper constructed an adaptive spline-based DSM fusion method guided by the geometric structure characteristics of the building and established an adaptive-meshing spline function to mine the structural information according to the complexity of the building roof structure so as to guide the effective fusion of multiple DSMs and output the building DSM with high accuracy while maintaining the target structural characteristics.

First, the building segmentation method [28] is used to determine the building target in the image, and then the elevation value of the building target in the n groups of DSMs is determined according to the results of the stereo solution. For n DSMs, the points of the building area in the i^th DSM are represented as a set $S_i|\begin{array}{l}n\\ i=1\end{array}$, and the j^th member in $S_i|\begin{array}{l}n\\ i=1\end{array}$ is a vector $\left[p_x^{\left(i,j\right)},p_y^{\left(i,j\right)},p_h^{\left(i,j\right)}\right]$, which represents the horizontal coordinate, vertical coordinate, and elevation value of the j^th point, respectively. The information of each point in the DSM to be output is represented by a set $S$, whose j^th member is a vector $\left[p_x^j,p_y^j\right]$, representing the horizontal and vertical coordinates of the point, and whose final result requires the output of the corresponding $p_h^j$. Then, the elevation of any point of the DSM of the fused building is obtained.

The elevation values of the fusion results of n sets of DSMs need to satisfy two points, one is that the difference between the elevation of the corresponding approximate homonymous points before fusion is small enough, and the other is that they are as consistent as possible with the potential roof geometry of the building based on the spline function. To address the first point, the confidence level of each point in the DSM should be determined first. Similar to the method of confidence definition in the fusion method of surface DSM, if the stereo solutions of the near homonymous points of different DSMs are correct, they should have similar elevations. According to this assumption, this paper used the similarity of elevation of the near homonymous points as the confidence degree of the corresponding points. For the second point, the mesh partitioning of the spline function needs to be defined. To fit buildings with different complex structures, the strategy of iterative adaptive meshing was utilized to construct the proposed DSM fusion algorithm. Considering that most of the buildings are flat and pitched roofs, the grid is initialized as a low-scale grid. If the elevation deviation of the obtained fusion result is less than a given threshold, the current grid division is valid; otherwise, the grid scale is increased, and recursive operations are performed. Based on the above two objectives, the optimization problem of n DSMs fusion can be constructed as shown in Equation (5).

$$\begin{array}{l}\underset{P_h,{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1}{\min }{F}_{fusion}\left(P_h,{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right)=\\ {\displaystyle \sum _i{\displaystyle \sum _j{c}_j{\Vert p_h^j-p_h^{\left(i,Nea{r}_i\left(j\right)\right)}\Vert }_2^2}}+{\delta }_1\times {\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\Vert \left[a_0^j,a_1^j,b_1^j,-1\right]\times {\left[1,p_x^i,p_y^i,p_h^i\right]}^T\Vert }_2^2}}\\ s.t.\left[a_0^j,a_1^j,b_1^j\right]\times {\left[1,n_x^i,n_y^i\right]}^T=\left[a_0^k,a_1^k,b_1^k\right]\times {\left[1,n_x^i,n_y^i\right]}^T\forall j\in node\left(i\right),k\in node\left(i\right)\\ \left[a_1^j,b_1^j\right]\times {\left[l_x^i,l_y^i\right]}^T=\left[a_1^k,b_1^k\right]\times {\left[l_x^i,l_y^i\right]}^T\forall j\in line\left(i\right),k\in line\left(i\right)\end{array}$$

(5)

where $\left[a_0^j,a_1^j,b_1^j\right]$, $\left[1,n_x^i,n_y^i\right]$, $\left[l_x^i,l_y^i\right]$, $node\left(i\right)$ and $line\left(i\right)$ are the parameters of the j^th spline function, the gridline intersection coordinates, the gridline slope parameter, the spline label adjacent to node i, and the spline label adjacent to the gridline, respectively. The left half of the objective function represents the elevation residuals of the DSM before and after minimizing the fusion, and the right half represents the elevation residuals of the DSM and the spline function after minimizing the fusion. Constraints are used to control the continuity of the spline function at the grid lines. The above equation is a nonconvex optimization problem, and in order to solve the problem, the idea of alternating iterations is used to fix the variable ${\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T$ when solving for the variable $P_h$. At this time, the constraints do not contain variables, and the optimization problem is an unconstrained convex optimization problem, which can be obtained from the equation about $P_h$ by setting the gradient to 0 by the objective function, as shown in Equation (6).

$$\frac{1}{N}{\displaystyle \sum _{i}diag\left(C\right)\left[P_h-P_h^i\right]}-{\delta }_1\times {\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}\left[a_0^j,a_1^j,b_1^j,-1\right]\times {\left[1,p_x^i,p_y^i,p_h^i\right]}^T}}=0$$

(6)

Therefore, the optimal solution $P_h$ can be expressed as Equation (7).

$$P_h^{\ast }={\left(\frac{1}{N}{\displaystyle \sum _{i}diag\left(C\right)}-diag\left(-1_N\right)\right)}^{-1}\left({\delta }_1\times {\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}\left[a_0^j,a_1^j,b_1^j\right]\times {\left[1,p_x^i,p_y^i\right]}^T}}+\frac{1}{N}{\displaystyle \sum _{i}diag\left(C\right)\left[P_h^i\right]}\right)$$

(7)

Subsequently, the obtained $P_h^{\ast }$ is fixed and solved for the variable ${\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T$. At this point, the optimization problem is an equation-constrained convex optimization problem, and the original problem is first transformed into a Lagrangian function to remove the equation constraint, as shown in Equation (8).

$$\begin{array}{l}L\left({\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1,\gamma \right)=\\ {\displaystyle \sum _i{\displaystyle \sum _j{c}_j{\Vert p_h^j-p_h^{\left(i,Nea{r}_i\left(j\right)\right)}\Vert }_2^2}}+{\delta }_1\times {\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\Vert \left[a_0^j,a_1^j,b_1^j,-1\right]\times {\left[1,p_x^i,p_y^i,p_h^i\right]}^T\Vert }_2^2}}+\gamma A{\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T\end{array}$$

(8)

where $A{\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T=0$ is the collapsed equation constraint, $\gamma$ is the dual variable. Further, the Lagrangian dual function is shown in Equation (9).

$$\begin{array}{l}G\left(\gamma \right)=\underset{P_h}{\min }{\displaystyle \sum _i{\displaystyle \sum _j{c}_j{\Vert p_h^j-p_h^{\left(i,Nea{r}_i\left(j\right)\right)}\Vert }_2^2}}+\\ {\delta }_1\times {\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\Vert \left[a_0^j,a_1^j,b_1^j,-1\right]\times {\left[1,p_x^i,p_y^i,p_h^i\right]}^T\Vert }_2^2}}+\gamma A{\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T\end{array}$$

(9)

According to the slater condition, when the objective function is quadratic programming and the constraints are linear constraints, the dual gap is 0, and the optimal solution of the original problem coincides with the optimal solution of the dual problem. Since the KKT condition is a sufficient condition for obtaining the optimal solution of the Lagrangian dual problem, the optimal solution of the dual problem can be solved through the KKT condition, as shown in Equation (10).

$$\begin{array}{l}\frac{\partial L}{\partial \left[{\mathit{a}}_0^{\ast },{\mathit{a}}_1^{\ast },{\mathit{b}}_1^{\ast }\right]}=0\Rightarrow \\ {\delta }_1\times {{\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\left[1,p_x^i,p_y^i\right]}^T\times \left[1,p_x^i,p_y^i,p_h^i\right]\times \left[a_0^{j*},a_1^{j*},b_1^{j*},-1\right]}}}^T+{\gamma }^{T}A=0\\ A{\left[{\mathit{a}}_0^{\ast },{\mathit{a}}_1^{\ast },{\mathit{b}}_1^{\ast }\right]}^T=0\end{array}$$

(10)

To solve for the optimal $\left[{\mathit{a}}_0^{\ast },{\mathit{a}}_1^{\ast },{\mathit{b}}_1^{\ast }\right]$, the above equation is replaced by $\left[{\mathit{a}}_0^{\ast },{\mathit{a}}_1^{\ast },{\mathit{b}}_1^{\ast }\right]=\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]+\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1\right]$ and rewritten in matrix equation form (see Equation (11)) to further solve for $\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1\right]$.

$$\begin{array}{l}\left[\begin{array}{l}H,A^T\\ A,0\end{array}\right]{\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1,\gamma \right]}^T=\left[\begin{array}{l}-G\\ 0\end{array}\right]\\ G={\delta }_1\times {{\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\left[1,p_x^i,p_y^i\right]}^T\times \left[1,p_x^i,p_y^i,p_h^i\right]\times \left[a_0^j,a_1^j,b_1^j,-1\right]}}}^T\\ H={\delta }_1\times {\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\left[1,p_x^i,p_y^i\right]}^T\times \left[1,p_x^i,p_y^i\right]}}\end{array}$$

(11)

Solve this matrix equation to obtain $\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1\right]$, and update the optimization variables according to $\left[{\mathit{a}}_0^{\ast },{\mathit{a}}_1^{\ast },{\mathit{b}}_1^{\ast }\right]=\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]+\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1\right]$, alternately optimize $P_h$ and ${\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T$, until $\left|F_{fusion}\left(P_h^n,{\mathit{a}}_0^n,{\mathit{a}}_1^n,{\mathit{b}}_1^n\right)-F_{fusion}\left(P_h^{n-1},{\mathit{a}}_0^{n-1},{\mathit{a}}_1^{n-1},{\mathit{b}}_1^{n-1}\right)\right|\le T_{threshold}$ is satisfied or the number of iterations exceeds a given threshold. When the iterations converge, the average elevation deviation between the fused DSM and the spline function is compared, that is, $\frac{1}{N_{DSM}}{\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\Vert \left[a_0^j,a_1^j,b_1^j,-1\right]\times {\left[1,p_x^i,p_y^i,p_h^i\right]}^T\Vert }_2^2}}$, where $N_{DSM}$ is the number of points in the fused DSM. If the average elevation deviation is less than the threshold, the current scale of meshing can effectively describe the building roof structure; otherwise, the scale of meshing is increased until the elevation deviation is less than the threshold or the number of iterations exceeds the limit. The DSM fusion method based on adaptive spline guided by building target structure characteristics can be summarized in Algorithm 1.

Algorithm 1. DSM fusion method based on adaptive spline guided by building target structure characteristics.

Input: n sets of DSMs $S_i|\begin{array}{l}n\\ i=1\end{array}$   Initialization: Initialize the mesh partitioning of the spline function, calculate the n sets of DSM weights $C$ and the iteration counter $n_1=0,n_1=0$, initialize $P_h^n$ and ${\left[{\mathit{a}}_0^n,{\mathit{a}}_1^n,{\mathit{b}}_1^n\right]}^T$.   Step 1: $n_1=n_1+1$   Step 2: $n_2=n_2+1$   Step 3: Use alternate iteration method to optimize variables, fix ${\left[{\mathit{a}}_0^n,{\mathit{a}}_1^n,{\mathit{b}}_1^n\right]}^T$ and $P_h^n$ in turn, and use Equations (3) and (7) to calculate $P_h^{\ast }$ and $\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1\right]$.   Step 4: Let $P_h^n=P_h^{\ast }$, ${\left[{\mathit{a}}_0^n,{\mathit{a}}_1^n,{\mathit{b}}_1^n\right]}^T={\left[{\mathit{a}}_0,{\mathit{a}}_1,{\mathit{b}}_1\right]}^T+\left[\Delta {\mathit{a}}_0,\Delta {\mathit{a}}_1,\Delta {\mathit{b}}_1\right]$.   Step 5: If $\left|F_{fusion}\left(P_h^n,{\mathit{a}}_0^n,{\mathit{a}}_1^n,{\mathit{b}}_1^n\right)-F_{fusion}\left(P_h^{n-1},{\mathit{a}}_0^{n-1},{\mathit{a}}_1^{n-1},{\mathit{b}}_1^{n-1}\right)\right|\le T_{threshold}$ or $n_2\ge n_{th}^2$ then execute Step 6; otherwise, skip to Step 2.   Step 6: If $\frac{1}{N_{DSM}}{\displaystyle \sum _j{\displaystyle \sum _{\left[p_x^i,p_y^i\right]\in S_j}{\Vert \left[a_0^j,a_1^j,b_1^j,-1\right]\times {\left[1,p_x^i,p_y^i,p_h^i\right]}^T\Vert }_2^2}}\le h_{th}$ or $n_1\ge n_{th}^1$ then the algorithm is executed; otherwise, increase the mesh density and skip to Step 1.   Output: DSM after fusion $S=\left\{P^{n_2}|1\le i\le N_{DSM}\right\}$.

3. Results

3.1. Experimental Results and Analysis of Simulated DSM Data

In order to verify the performance of the fusion algorithm under different roof shape building conditions, firstly, the target DSM with different roof shapes was generated as the truth map using computer simulation, and then Gaussian noise with different variances was generated and superimposed on the truth map DSM to generate two sets of observations, called DSM 1 and DSM 2. The two sets of noise-doped DSM 1 and DSM 2 were fused by executing the proposed DSM fusion algorithm, and the fusion results are shown in Figure 4, Figure 5 and Figure 6.

It can be seen from the figure that after the noise is superimposed, the two sets of observed inputs DSM1 and DSM2 produce different degrees of difference from the truth map. When the noise variance is greater than 0.5, the roof structure of the building has undergone major deformation. By executing the proposed adaptive spline-based target characteristic-guided DSM fusion algorithm, the potential roof shapes in the DSMs to be fused are explored to the maximum extent to guide the generation of DSMs with high confidence. It can be seen from Figure 5d,g,j that the fused DSMs have the same roof geometry as the truth map, especially when the noise variance is 1. The proposed method still obtains excellent fusion results, and the comparison of elevation errors before and after fusion is shown in

Table 1. Comparison of elevation errors before and after fusion of flat-roofed buildings.

	DSM 1	DSM 2	DSM after Fusion
Noise variance 0.5	0.3934	0.3976	0.0128
Noise variance 1	0.8120	0.7076	0.0135

,

Table 2. Comparison of elevation errors before and after fusion of pitched-roofed buildings.

	DSM 1	DSM 2	DSM after Fusion
Noise variance 0.1	0.0810	0.0801	0.0762
Noise variance 0.5	0.3942	0.4076	0.1266
Noise variance 1	0.8226	0.7983	0.1268

and

Table 3. Comparison of elevation errors before and after fusion of four-slope-roofed buildings.

	DSM 1	DSM 2	DSM after Fusion
Noise variance 0.5	0.3858	0.4095	0.0203
Noise variance 1	0.7961	0.7609	0.0320

.

The elevation error calculation is shown in Equation (12).

$$\sqrt{\frac{1}{m}{\displaystyle \sum _{i=1}^m{(h_i-{\widehat{h}}_i)}^2}}$$

(12)

where $h_i$ represents the true elevation value of the ground point, ${\widehat{h}}_i$ represents the generated elevation value, and m represents the total number of ground points.

Through Table 1, Table 2 and Table 3, it can be found that the proposed DSM fusion algorithm can still significantly improve the DSM accuracy of the target under different roof shape building conditions, even in the case of large noise variance.

3.2. IKONOS Data Experiment Results

In this paper, two different stereo matching algorithms are used to generate two sets of DSMs for IKONOS data: Method 1 is a stereo matching algorithm based on the modified Census transform [29], Method 2 is a robust stereo matching algorithm using adaptive random walk with restart [30], and then the proposed DSM fusion algorithm is executed to obtain the final DSM. Combined with the truth map to compare and analyze the errors, the experimental results are shown in Figure 7, and the RMSEs before and after DSM fusion are shown in

Table 4. Comparison of RMSE before and after fusion.

Method	RMSE
Method 1	2.24
Method 2	2.87
Fusion	2.16

. The experimental results showed that the proposed method can effectively constrain and improve the DSM of the building, and the integrity of the overall construction of the target 3D model structure was significantly improved.

4. Discussion

In this paper, we studied the DSM fusion method for building targets. In order to obtain better DSM results for buildings, a novel building DSM fusion method based on adaptive splines and target characteristic guidance was proposed. The improved DSM was obtained by assigning different DSM weights to integrate the advantages of multiple DSMs. Based on the existing prior knowledge, the better the DSM effect, the higher the weight assigned; otherwise, it is treated as equal weight.

This paper used adaptive spline functions to fit the roof structure of different buildings. Since the spline function is a piecewise function, it can be divided into an appropriate interval to accurately describe the specific building roof structure. At the same time, through rigorous mathematical proofs, we concluded that the denser the grid, the more accurately the spline function can describe the complex roof shape of a building. Therefore, this paper constructed a set of adaptive-meshing low-order spline functions to describe the characteristics of the target structure. The adaptive spline method can be applied in many aspects [31,32,33].

The part of the local deformation of the building in the DSM generated by the inaccurate elevation will also be fitted by the spline function. The DSM obtained directly is not accurate enough, and many existing studies are based on the constraint DSM [1,3,5,27,34]. In order to maintain the building roof structure as much as possible and avoid the introduction of inaccurate elevation information in the DSM, a DSM fusion method based on adaptive spline guided by building geometric structure characteristics was constructed to obtain a more accurate DSM of the building by fusing multiple DSMs. The acquired DSMs can come from different sources, either from LIDAR or from stereo reconstruction [29,30]. We took the fusion of two DSMs as an example, used simulation to generate DSMs with different roof shapes, generated two sets of observation data by superimposing Gaussian noise with different variances in the generated DSMs, then performed fusion experiments using the proposed method. At the same time, we also used IKONOS data with different stereo matching algorithms to generate the corresponding DSMs and perform fusion. The greater the number of DSMs, the better the fusion effect will be. However, as the number of DSMs increases, the corresponding computational consumption will also increase. In order to balance the number of DSMs with computational consumption, the number of DSMs is not as large as possible.

In summary, the proposed method can effectively preserve the edges and the integrity of the DSM of buildings. The experimental results show that the proposed method is an effective and efficient DSM constraint method for buildings.

5. Conclusions

This paper proposed an adaptive spline and target characteristic-guided DSM fusion update method for buildings and conducted DSM fusion experiments using typical urban area images of different scenes. The experimental results showed that the proposed method can effectively constrain and improve the DSM of buildings, and the integrity of the overall construction of the target 3D model structure was significantly improved, which is an effective and efficient DSM constraint method for buildings.