Comparative Analysis of Digital Elevation Model Generation Methods Based on Sparse Modeling

Abstract

With the spread of aerial laser bathymetry (ALB), seafloor topographies are being measured more frequently. Nevertheless, data deficiencies occur owing to seawater conditions and other factors. Conventional interpolation methods generally need to produce digital elevation models (DEMs) with sufficient accuracy. If the topographic features are considered as a basis, the DEM should be reproducible based on a combination of such features. The purpose of this study is to develop new DEM generation methods based on sparse modeling. Based on a review of the definitions of sparsity, we developed DEM generation methods based on a discrete cosine transform (DCT), DCT with elastic net, K-singular value decomposition (K-SVD), Fourier regularization, wavelet regularization, and total variation (TV) minimization, and conducted a comparative analysis. The developed methods were applied to artificially deficient DEM and ALB data, and their accuracy was evaluated. Thus, as a conclusion, we can confirm that the K-SVD method is appropriate when the percentage of deficiencies is low, and that the TV minimization method is appropriate when the percentage of deficiencies is high. Based on these results, we also developed a method integrating both methods and achieved an RMSE of 0.128 m.

1. Introduction

Digital elevation models (DEMs) comprise the primary data structure for representing detailed topographic information and are used in various ways [1]. Aerial laser surveying and aerial laser bathymetry (ALB) have been widely used to create DEMs. For example, ALB data can be used for various purposes [2], such as in sediment volume surveys of dams and rivers and seafloor surveys for offshore construction. However, in these measurements, cases occur where measurement data are missing owing to topographic surface conditions, or obstructions [3]. In ALB, the performance of bathymetry can be degraded owing to various factors such as the water quality and surface and bottom conditions. Moreover, the density of the laser measurement data in a water area is approximately 1/10 of that in a land area. Generally, interpolation is used to generate a DEM from deficient data. Currently, interpolation methods such as the nearest neighbor method, Kriging, and inverse distance weighting are used to generate DEMs [4]. However, their accuracy can degrade in regions with sparse data. When interpolation with a guaranteed accuracy cannot be performed, re-measurement can become necessary owing to the low quality of the interpolation, creating a severe problem. Therefore, a high-accuracy method for interpolating DEMs is required.

It is worthwhile to consider the characteristics of the physical shape of the topography. Topographic quantities such as slopes, curvatures, and ground clearances have been commonly used to describe topographic features. These topographic features have been used to analyze the geomorphology of a target area by introducing knowledge of topographic formation mechanisms. In addition, many topographic classifications and representations have been created based on topographic features [5]. These results indicate that even limited topographic features can effectively represent the topography. If the topographic features are used as a basis, a combination of such features should make it possible to reproduce the topography. However, in the context of topographic reconstruction based on limited features, it remains unclear which bases are generally effective. In addition, from the perspective of generalizability, it is desirable to use as few bases as possible. Sparse modeling is a method that satisfies these requirements. It can extract important bases by assuming sparsity even when the number of data is limited. Owing to these characteristics, significant results have been achieved in image restoration. For example, the K-singular value decomposition (K-SVD) method has attracted attention due to its assumption of sparsity in the image space [6]. K-SVD uses a set of many bases (referred to as a dictionary) to represent an image. It assumes that only a small number of bases from this dictionary can be used to represent the image. The details of K-SVD are described in Section 3.2.2. Lustig et al. used an interpolation method that assumes sparsity in the frequency space of the image [7]. Here, the assumption is that an image can be represented using a small number of various frequency patterns.

As described above, various definitions of sparsity exist in image restoration. Similarly, although various types of sparsity can be assumed for topography, understanding which type of sparsity is effective for the topography requires further clarification. Therefore, it is necessary to construct several generation methods based on sparsity for topography and to conduct a comparative analysis. Moreover, the assumption of sparsity in topography provides a perspective not found in conventional interpolation methods, as it interpolates by considering the spatial correlations of topographic features.

This study aims to develop such new DEM generation methods. They should enable efficient DEM generation and provide high-resolution information. In order to achieve this aim, this study constructs multiple DEM generation methods by considering various types of sparsity in sparse modeling and performs a comparative analysis. Specifically, the definitions of sparsity in topography are compared, and multiple topography interpolation methods corresponding to each sparsity definition are considered. Then, each method is applied to artificially deficient data from an existing DEM to determine the topography interpolation performance based on the percentage of deficient data. Then, the applicability of each method is verified on actual data obtained from ALB measurements. Finally, we discuss the results and summarize the characteristics of each method.

2. Related Work

DEM generation can be classified into methods for generating DEMs from nodal data arranged in a grid of points, and methods for generating DEMs from randomly arranged point data. In the former, many DEM generation methods assume self-similarity, namely, the existence of topographic similarity in the shape of the topography at a given point in the surrounding area. One popular type uses filters (e.g., median or nonlinear filters) while considering the spatial similarity [8]. However, the applicability of these methods has only been verified in cases where approximately 2% of the total data, at most, are missing. Other methods include spline interpolation [9], inverse filtering [10], and singular value decomposition [11]. Additionally, “absolutely minimizing Lipschitz extension” has been proposed; this is an interpolation method for addressing missing regions or elevation discontinuity boundaries caused by sea levels or inland water surfaces in the DEM [12]. However, none of these methods have been validated for cases where the missing data are located in the entire area of interest.

The methods for generating DEMs from randomly arranged point data include the use of weighted average and weighted first-order interpolation methods, in which the values around each point are approximated by a first-order plane [13]. Another generation method determines the surface while aiming to minimize the second-order variance representing the smoothness of the surface [14]. Nevertheless, owing to their deficiencies, these approaches have limited applicability to cases with differences in the point density.

As mentioned above, this study focuses on sparse modeling [15]. Sparse modeling was separately developed for sparse coding [16] and compressed sensing [17]. Sparse coding was developed mainly for signal processing, such as for images and audio. It is used for various applications such as noise reduction, interpolation, and classification [18]. The concept of sparse coding is that “a natural signal is approximated by using a large set of bases (a dictionary) and combining only a small number of bases from this dictionary”. In other words, sparse coding aims to find a sparse representation of the observed data using a given dictionary. Various algorithms have been proposed for this problem. One example is orthogonal matching pursuit (OMP) [19,20]. This algorithm is a greedy method that individually adds the best basis candidates rather than selecting the best basis from among all combinations [21]. The details of OMP are described in Section 3.2.1. “Basis Pursuit” is another example [22]. This method is attributed to the relaxation problem of minimizing the L₁ norm instead of the L₀ norm in the sparse representation. A method called the “Least Absolute Shrinkage and Selection Operator” adds a regularization term for the L₁ norm to the error function of the ordinary least squares method [23,24,25]. Overall, however, the accuracy of such sparse coding methods depends on how well a given dictionary fits the data.

Two methods are mainly used for determining dictionaries. The first method uses an existing dictionary. Examples include the discrete cosine transform (DCT), wavelet transform, curvelet transform, and bandelet transform. Starck et al. and Rubinstein et al. described the above dictionaries in detail, and we refer the reader to those studies for the details [26,27]. The second method involves learning a dictionary. In one example, two steps are alternated to perform the dictionary learning. First, sparse coding uses a given dictionary to identify sparse representations of the observed data. Second, the dictionary is updated by fixing the values of the sparse representation coefficients and positions of nonzero elements. Olshausen et al. proposed an algorithm using these two steps for dictionary learning [28,29]. Engan et al. improved on the high computational cost of that algorithm and proposed a more efficient algorithm called the “Method of Optimal Directions” [30]. This method requires the computation of inverse matrices to update the dictionary, making it relatively computationally expensive. Aharon et al. proposed the K-SVD method [6]. In this method, bases in the dictionary are updated one by one instead of updating the dictionary all at once [31]. The update is approximated by using SVD.

Compressed sensing, conversely, is a framework for reconstructing anonymous data from a small number of incomplete observations [32,33]. It assumes that the data are sparse in the original space or in a space after some transformation process. Compressed sensing uses this sparsity to enable reconstruction from a small number of samples (i.e., less than that from the sampling frequency). Compressed sensing is also based on random sampling. For example, consider the reconstruction of magnetic resonance images in a k-space. If the observed data are equally sampled in the k-space, the reconstruction results in a folded image called a folded artifact [34]. As the original signal and its duplicate are indistinguishable, it is impossible to reconstruct the original image from the observation data. When random sampling is used, the resulting artifacts appear as Gaussian noise. The incoherence of this noise makes compressed sensing possible. Compressed sensing was first introduced by Donoho et al. [35]. Since then, compressed sensing has been utilized in various fields and domains, such as in medicine and astronomy [7,36].

3. Methods

This section first summarizes DEM generation approaches using sparse modeling. Typical sparsity approaches include sparse coding, compressed sensing, and total variation (TV). We developed DEM generation methods according to the definitions of sparsity.

3.1. Definitions of Sparsity and Comparative Analysis Methods

Sparse coding uses a set of many bases (a dictionary) to represent the topography. However, only a few bases from this dictionary can be used to represent the topography. Here, the approach focuses on a local topography. When the coefficients are computed using the observed topography and dictionary (i.e., sparse coding), the sparsity of the coefficients is assumed. Sparsity is assumed such that any localized 2D plane cut from the topography can be represented by a linear sum of a small number of the many elements/bases comprising the topography (Figure 1). This approach focuses on spatial patterns and is based on generating the topography from a small number of spatial patterns.

Compressed sensing assumes that the topography can be represented using a small number of patterns among various frequency patterns. It assumes sparsity in the frequency space of the topography, that is, only some coefficients are significant, while others are zero or very small. It assumes that a section of the topography can be considered as a wave and expressed as a linear sum of a small number of patterns among the many elements and underlying waves (frequency patterns) comprising the topographic section (Figure 2). The topography is generated from a small number of frequency patterns after the entire terrain is transformed into the frequency space.

The TV approach considers that the topography is smooth and that the differences between adjacent points are often tiny, with few areas of significant change. Sparsity is assumed for such topographic variations. This approach is based on the sparsity of the few sizeable topographic change points.

Based on Section 2 and the above definitions, the following six methods were comparatively analyzed in this study.

• DCT method: This is a sparse coding framework that uses a redundant DCT dictionary (an existing general-purpose dictionary) instead of learning a dictionary. The DCT transforms a discrete sequence of signals into a combination of cosines with various frequencies and amplitudes, a technique derived from the discrete Fourier transform. The sparse coding is performed by using the dictionary to reproduce the topography.
• DCT with elastic net method: This is another sparse coding framework. It also uses a redundant DCT dictionary (an existing general-purpose dictionary) instead of learning a dictionary. The sparse coding is performed using the dictionary to reproduce the topography. The sparsity assumption is relaxed during the sparse coding, and elastic net regression is used. This allows for expanding the range of the local topography.
• K-SVD method: another sparse coding framework, the K-SVD method learns a dictionary from an observed deficient topography and performs sparse coding using the dictionary to reproduce the topography.
• Fourier regularization method: this compressed sensing framework assumes the sparsity of the topography in the Fourier space and reconstructs the topography.
• Wavelet regularization method: this compressed sensing framework assumes the sparsity of the topography in the wavelet space and reconstructs the topography.
• TV minimization method: this method reproduces the topography by minimizing the sum of the differences (gradients) between adjacent points (i.e., the TV).

3.2. Digital Elevation Model Generation Method Using Sparse Coding

As noted above, two types of methods are based on the sparse coding framework: those that use existing dictionaries and those that learn dictionaries. In both types, the DEM is represented as a linear sum of a dictionary A and its sparse coefficients x.

3.2.1. Using Existing Dictionaries

First, DEM data with defects are cut into n × n patches to obtain N patches as patch data. This is based on the assumption that in the n × n patch, observation $y_i$ can be expressed as a linear sum of a dictionary A and its sparse coefficients $x_i$. $x_i$ is obtained from the known $y_i$ and A by solving the optimization problem in Equation (1) (Figure 3).

$${\widehat{x}}_i=\underset{x_i}{\arg \min }\frac{1}{2}\Vert y_i-MA{x}_i\Vert {}_2^2,s.t.x_{i}issparse$$

(1)

In the above, M is a mask matrix. The algorithm for solving the optimization problem uses OMP or elastic net regression, as described below. As there are deficiencies in $y_i$, the computation process is performed only on the portion without deficiencies. Specifically, the mask matrix M is introduced to align the number of rows in each matrix. The dictionary is a redundant DCT dictionary, that is, a standard dictionary used in image compression such as JPEG (Figure 4). From the obtained sparse coefficients $x_i$ and dictionary A, the topography of the corresponding patch can be reproduced and interpolated as ${\widehat{y}}_i$ using Equation (2).

$${\widehat{y}}_i=A{x}_i$$

(2)

The value of each mesh in the DEM is obtained by averaging the values of the corresponding cells in all patches containing that mesh.

In OMP, Equation (1) is regarded as the following optimization problem; the solution is obtained using a greedy method.

$${\widehat{x}}_i=\underset{x_i}{\arg \min }\frac{1}{2}\Vert y_i-MA{x}_i\Vert {}_2^2,s.t.\Vert x_i\Vert {}_0\le k_0$$

(3)

where k₀ is a constant that defines the sparsity. First, x₀ = 0 and k = 0 are set as initial values. Initial values are also set for the initial residuals r⁰ = $y_i$, sparsity k₀, and initial solution support S⁰ = $\varphi$. Next, k is incremented by one, and the basis vector a_j that maximizes the inner product $\left|a_j·r^k\right|$ with the residual $r^k$ from the dictionary A is added to the support $S^k$. The tentative solution $x^k$ is then updated by solving Equation (4), and the residual $r^k$ is updated based on Equation (5). The above process is repeated k₀ times.

$$x^k=\underset{x_i}{\arg \min }\frac{1}{2}\Vert y_i-A{x}_i\Vert {}_2^2,s.t.Support\left\{x_i\right\}=S^k$$

(4)

$$r^k=y-A{x}^k$$

(5)

Elastic net regression, in contrast, is considered as a linear regression problem that uses the L₁ and L₂ norms as regularization terms for Equation (1) [37,38].

$${\widehat{x}}_i=\underset{x_i}{\arg \min }\frac{1}{2}\Vert y_i-MA{x}_i\Vert {}_2^2+{\lambda }_1\Vert x_i\Vert {}_1+{\lambda }_2\Vert x_i\Vert {}_2$$

(6)

where λ₁ and λ₂ are the weights of the L₁ and L₂ regularization terms, respectively. In elastic net regression, as with L₁ regularization, it is easy to obtain sparse solutions. In L₁ regularization, for example, only one solution is likely to be chosen when dealing with two variables. Elastic net regression, in contrast, makes it easy to select both variables simultaneously owing to the effect of the L₂ regularization term. In other words, the sparsity constraint is slightly relaxed. As discussed in more detail below, the patch size must be increased when considerable deficiencies exist. If the sparsity restrictions are excessively strict, the estimation may not be possible. Therefore, in this case, sparse coding is performed using elastic net regression as in Equation (6).

3.2.2. Using Dictionary Learning

In this study, we focus on the K-SVD method [6] as a dictionary-learning method. As in the previous section, we describe the K-SVD method in the context of DEM data with deficiencies.

First, the defective DEM data are cut into n × n patches to obtain N patches of data. The initial values of k = 0 and the sparsity k₀ are set, and a redundant DCT dictionary is set as the initial dictionary A₀. As the defective DEM data for the dictionary learning, we randomly sample from the N patches and use them to create $y_k$. The sparse representations $x_k$ are estimated from the data $y_k$ and dictionary $A_k$ by solving Equation (7) using OMP described above. The sum of products of each basis and the corresponding sparse representation are expressed as shown in Equation (8). Figure 5 shows an image of the sum of products. The bases of the dictionary are updated one by one. As an example, let $a_i$ be the basis to be updated. The sum of products of this basis and the corresponding sparse representation can be expressed as shown in Equation (9) using the data $y_k$, basis not to be updated $a_{j\ne i}$, and corresponding sparse representation $x_{\mathit{j}\ne \mathit{i}}^{\mathit{T}}$. The updated $a_{\mathit{i}}{x}_{\mathit{i}}^{\mathit{T}}$ can be obtained by performing a low-rank approximation of $a_{\mathit{i}}{x}_{\mathit{i}}^{\mathit{T}}$ using singular value decomposition. Then, the sequential updating of all of the bases updates the dictionary.

$$\widehat{x}=\underset{x_k}{\arg \min }\frac{1}{2}\Vert y_k-M{A}_k{x}_k\Vert {}_2^2,s.t.\Vert x_k\Vert {}_0\le k_0$$

(7)

$$A_{\mathit{k}}{x}_{\mathit{k}}\approx {\displaystyle \sum }_{\mathit{j}}{a}_{\mathit{j}}{x}_{\mathit{j}}^{\mathit{T}}$$

(8)

$$a_{\mathit{j}}{x}_{\mathit{j}}^{\mathit{T}}\approx y_{\mathit{k}}-{\displaystyle \sum }_{\mathit{j}\ne \mathit{i}}{a}_{\mathit{j}}{x}_{\mathit{j}}^{\mathit{T}}$$

(9)

3.3. Digital Elevation Model Generation Method Using Compressed Sensing

The Fourier and wavelet transforms [39] map the original signal to a frequency space. In compressed sensing, the data’s Fourier transform coefficients or wavelet transform coefficients are used as regularization terms to set up the optimization problem, as shown in Equation (10).

$$\widehat{x}=\underset{x}{\arg \min }\frac{1}{2}\Vert y-\mathrm{Mx}{\Vert }_2^2+\lambda \Vert \mathsf{\Psi }x{\Vert }_1$$

(10)

where $\mathsf{\Psi }$ represents the sparse transformation of each transform, and λ is the weight of the regularization term. To solve this optimization problem, we consider the more general optimization problem for non-smooth convex functions as an asymptotic solution (Equation (11)).

$$\underset{x}{\min }F\left(x\right)\equiv f\left(x\right)+g\left(x\right)$$

(11)

where $f\left(x\right)$ is a continuous and differentiable convex function, and $g\left(x\right)$ is a continuous convex function that is not necessarily smooth. In this case, the solution x is obtained based on the proximity gradient method, as follows.

First, the solution’s initial value x₀ and step width ρ are set. Using Equation (12), the gradient method is updated when only the function $f\left(x\right)$ is used. Next, the solution is updated using the proximity map ${\mathrm{prox}}_{{\lambda }_g}\left(r\right)$ according to Equation (13). This proximity map is given by Equation (14). The solution $\widehat{x}$ is obtained by iteratively applying Equations (12) and (13) until convergence. This method utilizes the structure of the function $g\left(x\right)$ and is called the “Iterative Shrinkage Thresholding Algorithm”, which is especially useful for L₁ regularization.

$$r_{k+1}={\mathit{x}}_k-\rho \nabla f\left({\mathit{x}}_k\right)$$

(12)

$$x_{k+1}={\mathrm{prox}}_{{\lambda }_g}\left(r_{k+1}\right)$$

(13)

$${\mathrm{prox}}_{{\lambda }_g}\left(r\right)=\underset{x}{\arg \min }\frac{1}{2}\Vert x-r{\Vert }_2^2+\lambda g\left(x\right)$$

(14)

3.3.1. Compressed Sensing Using Fourier Transform

Assuming sparsity in the Fourier space in Equation (11) results in Equations (15) and (16).

$$f\left(x\right)=\frac{1}{2}\Vert y-\mathbf{Mx}{\Vert }_2^2$$

(15)

$$g\left(x\right)=\lambda \Vert \mathbf{F}\mathit{x}{\Vert }_1$$

(16)

where the function F is the function performing the Fourier transform. When the regularization term is expressed as shown in Equation (16), the proximity map becomes Equation (17).

$${\mathrm{prox}}_{{\lambda }_g}\left(r\right)=F^{T}soft\left(Fr,\lambda \right)$$

(17)

where the function $soft\left(y,\lambda \right)$ is called the “soft threshold function” and is defined according to Equation (18).

$$soft\left(y,\lambda \right)=\left\{\begin{array}{l}y-\lambda \left(\lambda <y\right)\\ 0\left(-\lambda \le y\le \lambda \right)\\ y+\lambda \left(y<-\lambda \right)\end{array}\right.$$

(18)

The solution can be obtained using the proximity gradient method as follows. First, the initial solution value x₀, step width α, and threshold value λ are set. Using Equation (19), the gradient method is used to update the solution. Next, the solution is updated by proximity mapping using Equation (20). The solution is updated only for defective meshes using Equation (21). In Equation (21), j refers to the j-th mesh, and y represents the defective DEM data. These steps are repeated until convergence to obtain the solution.

$$r_{k+1}={\mathit{x}}_k-\rho M^T\left(y-M{x}_k\right)$$

(19)

$$z_{k+1}=F^{T}soft\left(F{r}_{\mathit{k}+\mathbf{1}},\lambda \right)$$

(20)

$$x_{k+1}\left[j\right]=\left\{\begin{array}{c}{z}_{k+1}\left[j\right]\left(y\left[j\right]=NaN\right)\\ y\left[j\right]\left(y\left[j\right]\ne NaN\right)\end{array}\right.$$

(21)

3.3.2. Compressed Sensing Using Wavelet Transform

When assuming sparsity in the wavelet space, Equations (16) and (20) can be replaced by Equations (22) and (23), respectively. A similar algorithm can be used to obtain the solution.

$$g\left(x\right)=\lambda \Vert W\mathit{x}{\Vert }_1$$

(22)

$${\mathit{z}}_{k+1}={\mathit{W}}^{T}soft\left(\mathit{W}{\mathit{r}}_{\mathit{k}+\mathbf{1}},\lambda \right)$$

(23)

where the function W is the function performing the wavelet transform.

3.4. DEM Generation Method Using Total Variation Minimization

In this section, we first define the TV and explain the derivatives required for optimization. Then, a DEM generation method is presented based on the TV minimization criterion.

3.4.1. Total Variation

The TV norm is defined as the L₁ norm of the magnitude of the gradient in Equation (24) [40]. In the case of discrete data x like in a DEM, the TV norm can be defined as shown in Equation (25). Thus, the TV norm can be calculated by differentiating the DEM in the $x$ and $y$ directions, finding the gradients, and then summing them.

$$TV=\int \int \Vert \nabla f\left(x,y\right){\Vert }_{1}dxdy=\int \int \sqrt{{\left(\frac{\partial f}{\partial x}\right)}^2+{\left(\frac{\partial f}{\partial y}\right)}^2}dxdy$$

(24)

$$\Vert x_{TV}\Vert {}_1={\displaystyle \sum }_i{\displaystyle \sum }_j\left|\nabla x_{i,j}\right|={\displaystyle \sum }_i{\displaystyle \sum }_j\sqrt{{\left(x_{i,j}-x_{i-1,j}\right)}^2+{\left(x_{i,j}-x_{i,j-1}\right)}^2}$$

(25)

The derivative of the TV norm is expressed as shown in Equation (26). Equation (26) is not differentiable when the denominator is zero, but when a small quantity ε is introduced to approximate $\left|x\right|\approx \sqrt{x^{T}x+\epsilon }$, the derivative can be calculated as shown in Equation (27). For the actual calculations, the approximation in Equation (28) has been proposed [41]. Here, ε is a minimal value of approximately 10⁻⁸.

$$\nabla \Vert x_{TV}\Vert {}_1=\frac{\partial \Vert x_{TV}\Vert {}_1}{\partial x_{i,j}}$$

(26)

$$\frac{\partial \left|x\right|}{\partial x}\approx \frac{x}{\sqrt{x^{T}x+\epsilon }}$$

(27)

$$\begin{array}{ll}\frac{\partial \Vert x_{TV}\Vert {}_1}{\partial x_{i,j}}& \approx \frac{x_{i,j}-x_{i-1,j}}{\sqrt{{\left(x_{i,j}-x_{i-1,j}\right)}^2+{\left(x_{i-1,j+1}-x_{i-1,j}\right)}^2+\epsilon }}\\ & +\frac{x_{i,j}-x_{i,j-1}}{\sqrt{{\left(x_{i+1,j-1}-x_{i,j-1}\right)}^2+{\left(x_{i,j}-x_{i,j-1}\right)}^2+\epsilon }}\\ & -\frac{x_{i+1,j}-x_{i,j+1}-2x_{i,j}}{\sqrt{{\left(x_{i+1,j}-x_{i,j}\right)}^2+{\left(x_{i,j+1}-x_{i,j}\right)}^2+\epsilon }}\end{array}$$

(28)

3.4.2. Digital Elevation Model Generation Method Using Total Variation Minimization

Using the TV norm $\Vert x_{TV}\Vert {}_1$ from the previous section as the regularization term, we set up the optimization problem shown in Equation (29).

$$\widehat{x}=\underset{x}{\arg \min }\frac{1}{2}\Vert y-Mx{\Vert }_2^2+\lambda \Vert x_{TV}\Vert {}_1$$

(29)

For this optimization, we use the gradient method. First, the step width α is set. The derivative of the TV norm $\frac{\partial \Vert x_{TV}\Vert {}_1}{\partial x_{i,j}}$ is obtained using Equation (28). Only the solution for the defective mesh is updated according to Equations (30) and (31). j refers to the j-th mesh, and y represents the defective DEM data. These steps are repeated until convergence to obtain the solution.

$$z_{k+1}=x_k-\alpha \frac{\partial \Vert x_{TV}\Vert {}_1}{\partial x_{i,j}}$$

(30)

$$x_{k+1}\left[j\right]=\left\{\begin{array}{c}{z}_{k+1}\left[j\right]\left(y\left[j\right]=NaN\right)\\ y\left[j\right]\left(y\left[j\right]\ne NaN\right)\end{array}\right.$$

(31)

4. Results

As described in this section, we verified the performance of each method for different deficiencies. Then, the methods were applied to actual data from ALB measurements.

4.1. Application to Different Deficiency Ratios

The methods were applied to DEM data with artificial deficiencies. The original DEM (Mt. Unzen, Nagasaki, Japan) comprised 1 m mesh data (512 × 512 m) obtained from an aerial laser survey (Figure 6). Two areas with different topographic shapes were selected. Area 1 (upper area) has a constant slope throughout (elevation 560 to 470 m) and several valley lines. The elevation of the valley floor is approximately 450 m. Area 2 (lower area) is a cliff section, with elevations above and below the cliff of 660 and 510 m, respectively. We prepared two methods for selecting deficient meshes: one randomly selected a mesh, and the other selected a mesh by considering spatial correlations. The deficiencies were set as 5%, 10%, 20%, 25%, 30%, 40%, 50%, 60%, 70%, 75%, 80%, and 90% in the cases of random deficiencies. Figure 7 shows an example of the DEM data with deficiencies. The black areas in the figures indicate the deficiencies. The root mean square error (RMSE) was used to evaluate the accuracy of the DEM generation.

Figure 8 and Figure 9 show the application results. In the case of random deficiencies, the K-SVD method was effective when the percentage of deficiencies was under 50%. The DCT method was effective when the percentage was 50–75%, and the TV minimization method was effective when the percentage was 80% or more. Regarding the correlated deficiencies, the K-SVD and TV minimization methods were effective when the percentage of deficiencies was 75% or less. In comparison, the TV minimization method was effective when the percentage of deficiencies was 80% or more. The DCT and K-SVD methods had good accuracy when the percentage of deficiencies was low. Both methods used local features in the DEM, resulting in significant accuracy degradation when the percentage of deficiencies was high (80% and 60%, respectively). On the other hand, the DCT with elastic net method did not show the above significant accuracy degradation because it relaxed the constraint condition by using L₂ regularization. Compressed sensing, which deals with the frequency space of the entire DEM, did not show outstanding accuracy compared to the other methods. Comparison with the sparse coding results indicates that local features, rather than the entire DEM, are significant for terrain restoration. The TV method, on the other hand, considered local DEM gradients in the entire DEM, resulting in a higher overall accuracy. The effectiveness of each method did not differ depending on the topographic shape.

4.2. Application to Actual Aerial Laser Bathymetry Data

Next, each method was applied to actual ALB data. The data were measured in the shallow water area of a harbor (Port of Oniike, Kumamoto, Japan). As in the previous section, a 512 × 512 m section of the 1 m mesh data was cut out (Figure 10). The conditions of the seawater caused the deficiencies. The accuracy verification used bathymetric data at almost the same times. Cross-section 1 in the figure is a gentle valley with water depths ranging from 3 to 8 m. Significant deficiencies are observed, especially at the bottom of the valley. Cross-section 2 is a cliff section with water depths ranging from 4 to 9 m. ALB is capable of measuring at shallow depths, and this measurement confirms that there are many deficiencies, especially at water depths deeper than 7 m. The accuracy was verified by applying the DEM generation methods and calculating the RMSEs of the bathymetric data in cross-sections 1 and 2 in the figure. The nearest neighbor and Kriging methods were also compared as conventional methods.

Table 1. Root mean square errors (RMSEs) of applications to aerial laser bathymetry (ALB) data.

Methods	RMSE of Section 1 (m)	RMSE of Section 2 (m)
Discrete cosine transform (DCT)	0.463	0.257
DCT with elastic net	0.226	0.135
K-singular value decomposition (K-SVD)	10.837	0.296
Fourier regularization	0.353	0.292
Wavelet regularization	0.288	0.273
Total variation (TV) regularization	0.206	0.289
Nearest neighbor	0.526	0.296
Kriging	0.346	0.291

shows the application results. The effectiveness of each method did not differ depending on the topographic shape. This result appeared because both cross-sections had more deficiencies at their relatively flat bottoms.

5. Discussion

First, we discuss the results of Section 4.1 in the context of different deficiency ratios. The sparse coding methods without dictionary learning lead to ill-posed problems when the deficiency ratio is high; thus, the results could be more unstable. Elastic net can alleviate this instability because of its regularization effect. When using dictionary learning, the K-SVD method is effective when the deficiency ratio is low. This means that the K-SVD method can learn detailed topographies. In contrast, when the deficiency ratio is high, the dictionary cannot be appropriately learned owing to over-fitting, and the accuracy decreases. We also confirmed that the accuracy is improved when the dictionary is learned in a region with a small percentage of deficiencies and is shared across the entire range.

The differences between the sparse coding and compressed sensing methods are as follows. The sparse coding methods are applied to a local topography, whereas the compressed sensing methods are applied to the entire DEM in the frequency domain. Therefore, the compressed sensing method can generate a stable DEM even if localized areas have high deficiencies. The Fourier regularization method is more accurate than the wavelet regularization method for all deficiencies. From comparing the coefficients for the different regularization methods, it can be seen that the coefficients in the Fourier space are often smaller than those in the wavelet space. The Fourier regularization method is equally used for low, medium, and high frequencies. Therefore, the Fourier coefficients are also small, with Fourier coefficients below 0.001 accounting for 70%. The wavelet regularization method, on the other hand, uses mainly low-frequency components, with small wavelet coefficients accounting for only 40%. This is likely because the Fourier regularization represents global and local features. In addition, the TV minimization method was the most accurate and effective for all patterns, especially when the deficiencies were high.

Overall, the DEM generation results do not differ significantly depending on whether the defects are random or correlated. The accuracy of each method was also compared depending on the topography. The slope, Laplacian, and ground opening were used as quantitative indicators of the topography. No differences from the above trends were observed from any of these perspectives.

Next, we discuss the results from applying the methods to the ALB data as described in Section 4.2. The TV minimization method and DCT with elastic net method, which allow the patch size to be larger than the deficiencies, are effective in areas where the percentage of deficiencies is relatively high (line 1). The spatial extents of these deficiencies are also considered. In the case of small-scale deficiencies (line 2), the DCT with elastic net method is effective. This result arises from the effect of considering the smoothness in the L₂ regularization.

The results for the Fourier and wavelet regularization methods show relatively high RMSEs. All the data were summarized in the frequency space, which did not allow for capturing the local topography.

Based on the above discussion, we also developed a method integrating the K-SVD and TV minimization methods (TV-SVD method). The ALB data include areas with high and low percentages of deficiencies (low-density deficiencies). The practical application of the K-SVD method to the airborne laser DEM showed that the K-SVD method is effective for low deficiencies. In contrast, the TV minimization method is effective for high deficiencies.

The K-SVD results were not good because the estimation included areas with high deficiencies. The data in line 1 were divided into ranges with high and low deficiencies, and K-SVD was applied to the range with low deficiencies, resulting in a high accuracy (RMSE: 0.167 m). This suggests that dictionary learning works well only in areas with low deficiencies. It was also confirmed that a DEM could be generated accurately, even with shared dictionaries in different areas. This result shows the scalability of the dictionary and its ability to deal with many low deficiencies. Rather than learning dictionaries from data with high deficiencies, the DEM generation accuracy is higher when dictionaries are learned in different locations but in areas with low deficiencies.

In the case of large-scale deficiencies, the TV minimization method is effective because the dictionary learning does not work correctly. These results are consistent with previous results concerning the different deficiencies.

Accordingly, we applied the TV minimization method to the large-scale deficient parts and the K-SVD method to the small-scale deficient parts and integrated the two results. In order to construct a method according to the size of the deficiencies, it is necessary to distinguish the size. The large-scale deficient parts can be extracted using morphological operations [42]. Morphological operations include the closing operation, in which dilation is repeated a certain number of times, and then erosion is repeated the same number of times. The closing operation can be applied to a binary image in which deficient areas are set to 1 and other areas are set to 0 to extract regions. One of the characteristics of the closing operation is that it extracts large regions, while sparse regions are eliminated. The process makes it suitable for extracting large areas of deficiencies. It enables automatic identification of the applied method.

The proposed method was applied to the area around cross-section 1 in Figure 10. The morphological operation correctly identified the high-deficiency areas. Therefore, the K-SVD and TV minimization methods could be applied appropriately. The dictionary was learned in the same area. The RMSE of the result was 0.128 m. We confirm the improved accuracy over the results in Section 4.2. This result confirms the effectiveness of the proposed combination method. One of the results is shown in Figure 11. In Figure 11a, the uncolored areas are deficient data areas. Therefore, it looks as if a part of the figure is missing.

6. Conclusions

In this study, we compared and organized various types of sparsity in sparse modeling and constructed multiple DEM generation methods. In organizing the sparsity types, we focused on sparsity related to the basis of spatial topographic components, sparsity related to frequency patterns, and sparsity with few topographic change points. Based on the results, we developed DEM generation methods based on the DCT method, DCT with elastic net method, K-SVD method, Fourier regularization method, Wavelet regularization method, and TV minimization method. The developed methods were applied to artificially created DEMs with various deficiencies and to actual ALB data to verify the applicability of each method. In particular, we confirmed the effectiveness of the sparsity approach for spatial patterns and topographic changes. Specifically, we can confirm that the K-SVD method is appropriate when the percentage of deficiencies is low, and that the TV minimization method is appropriate when the percentage of deficiencies is high. Based on these results, we also developed a method integrating both methods and achieved an RMSE of 0.128 m. Furthermore, we clarified effective sparsity approaches for topography and effective DEM generation methods based on a discussion of the application results.

Future work includes extending the method to an integrated method for simultaneously considering multiple sparsities. In addition, the results from the multiple generation methods could be integrated into an ensemble. Such extended methods would be expected to stabilize the results and suppress the influence(s) of noise.