Sensitivity Analysis and Modeling for DEM Errors

Abstract

The Digital Elevation Model (DEM) can be created using airborne Light Detection And Ranging (LIDAR), Image or Synthetic-Aperture Radar (SAR) mapping techniques. The direct georeferencing of the DEM model is conducted using a GPS/inertial navigation system. The airborne mapping system datasets are processed to create a DEM model. To develop an accurate DEM model, all errors should be considered in the processing step. In this research, the errors associated with DEM models are investigated and modeled using Principal Component Analysis (PCA) and the least squares method. The sensitivity analysis of the DEM errors is investigated using PCA to define the significant GPS/inertial navigation data components that are strongly correlated with DEM errors. Then, the least squares method is employed to create a functional relationship between the DEM errors and the significant GPS/inertial navigation data components. The DEM model errors associated with airborne mapping system datasets are investigated in this research. The results show that the combined PCA analysis and least squares method can be used as a powerful tool to compensate the DEM error due to the GPS/inertial navigation data with about 27% in average for DEM errors produced by the direct georeferenced airborne mapping system.

1. Introduction

The Digital Elevation Model (DEM) can be created using airborne Light Detection And Ranging (LIDAR), Image or Synthetic-Aperture Radar (SAR) mapping techniques. The direct georeferencing of DEM model is conducted using GPS/inertial navigation system. The accuracy of the developed DEM model is strongly dependent on the mapping system and the georeferencing system grades. The selection of the mapping system and the georeferencing system grades is carried out in the planning stage based on the accuracy requirements of the required DEM model.

The previous literatures focused and heavily investigated the DEM model generation based on LIDAR, Images, and SAR data cleaning and filtering techniques such as Triangulated Irregular Network (TIN)-based filtering, slope-based filtering, mathematical morphological filtering, interpolation-based filtering, and Machine-learning-based filtering. The original TIN-based filtering was developed based on the classical progressive TIN densification (PTD) and was implemented effectively in the commercial software TerraScan [1, 2]. Then, the revised PTD was investigated and reduced the total errors by about 8% when compared with classical PTD method [3]. Afterwards, a Parameter-Free PTD (PFPTD) algorithm was developed and outperforms the classical and revised PTD methods [4]. The original slope-based filtering was derived based on height differences in the training dataset [5, 6]. Then, adaptive slope-based filtering algorithm was developed to improve the accuracy in urban applications when compared with the original slope-based filtering algorithm [7]. The original mathematical morphological filtering was proposed to filter LIDAR data [8]. Then, the progressive morphological filtering algorithm was developed to improve the original method by applying threshold condition based on the elevation differences and proposed increasing gradually the filtering window size [9]. Afterward, the spline iteration method was introduced to improve the morphological filtering algorithm [10]. All mathematical morphological filtering methods outputs are strongly dependent on adapting the filtering window size. The original interpolation-based filtering method was proposed to deal with the steep areas [11]. Then, the interpolation-based DEM generation method was developed where one of the Inverse Distance Weighted, Kriging, and Natural Neighbor (NN) can be employed for DEM generation [12]. The Natural Neighbor (NN) method was proven to provide most efficient results. Finally, Machine-learning-based filtering was investigated where this method depends on topographic characteristics of the areas under investigation [13]. The deep convolutional neural network (CNN) was proposed for develop accurate DEM model [14]. The Machine learning method was optimized using windowing method to improve the DEM model generation [15].

However, the above-mentioned methods are mainly dependent on the data cleaning and filtering techniques of the heights to develop the DEM model. These methods are associated with errors that are not considered in the DEM modeling and could be corelated with system navigation data. Therefore, this research investigated the sensitivity analysis modeling of DEM errors that are potentially correlated with navigation data to improve the overall accuracy of the DEM model.

2. Combined PCA and least squares method

The sensitivity and modeling of DEM errors are investigated in this research using PCA and least squares function modeling. The system navigation data (position, velocities, attitudes, accelerations and dopplers) are considered the inputs to the model and the DEM height error is considered the desired output of the model.

PCA is a numerical technique used to study multidimensional processes that can be used to (1) reduce the dimensionality of a dataset and (2) identify relationships between the underlying variables of the process. PCA is based on eigen or singular value analysis of the process correlation or covariance matrix. The goal of PCA is to determine the minimum number of eigenvectors that best describe the key features of the process correlation matrix. This results in a reduced-dimensionality model for the matrix which can be used for data analysis, reduction, and model synthesis. Singular value decomposition (SVD) is fundamental to PCA. More details on PCA and SVD can also be found in Jolliffe [16].

2.1 PCA analysis

Let X denote an m x n matrix. For convenience, we assume m ≥ n. The elements of the i^th row of X form the n-dimensional vector g_i. The elements of the j^th column of X form the m-dimensional vector a_j. The general singular value decomposition (gSVD) of X can be written as:

$$X=\mathit{U\Sigma }{V}^T$$E1

where U is an m x m matrix, Σ is an m x n matrix containing the singular values, and V^T is an n x n matrix. The columns of U are called the left singular vectors, {u_k}, and form an orthonormal basis for the range space, so that u_i·u_j = 1 for i = j, and u_i·u_j = 0 otherwise. The rows of V^T contain the elements of the right singular vectors, {v_k}, and form an orthonormal. U and V are orthonormal so that their inverses exist and are their transposes. The matrix Σ can be decomposed as:

$$\Sigma =\left[\begin{array}{c}S\\ 0\end{array}\right]$$E2

where S is an m x n diagonal matrix in which only the diagonal elements are non-zero, S = diag(s₁,…,s_n) where the diagonal elements are zero. If the rank of X is r, s_k > 0 for 1 ≤ k ≤ r, and s_i = 0 for (r + 1) ≤ k ≤ n. [For problems like the one we are interested in, noise generally ensures that r = n.] By convention, the ordering of the singular vectors is determined by high-to-low sorting of singular values, with the highest singular value in the upper left index of the S matrix. Note that for a square, symmetric matrix X, SVD is equivalent to eigenvalue decomposition. In PCA, the right singular vectors are frequently called the components. While the scaled left singular vectors {s_iu_k} are called the scores.

Note that U can be decomposed into two submatrices, an m x n matrix U_R and an m x m-n matrix U_N where U = [U_R U_N]. U_R defines the range space of U, while U_N defines the null space. Note that X = U_RSV^T so that XV = U_RS. This provides the reduced form of SVD often used in PCA. In practice this is the form generally used; hence, we often drop the R subscript on U. Figure 1 illustrates the various reduced-form matrices. Note that the right singular vectors span the space of the row vectors{g_i} and the left singular vectors span the space of the column vectors {a_j}.

Figure 1.

Illustration of reduced SVD matrices. The right singular vectors have often termed the components while the left singular values are the scaled scores.

Several relationships can be derived. The SVD equation for g_i is:

$${\mathit{g}}_{\mathit{i}}=\sum _{\mathit{k}=\mathbf{1}}^{\mathit{r}}{\mathit{u}}_{\mathit{ik}}{\mathit{s}}_{\mathit{k}}{\mathit{v}}_{\mathit{k}}$$E3

which is a linear combination of the right singular values {v_k}. The i^th row of U, g’_i, contains the coordinates of the i^th entry in the coordinate system (basis) of the scaled right singular values, s_kv_k. If r < n (or if we truncate the singular values to r = l), this computation requires fewer variables using g’_i rather than g_i, thus reducing the dimension of the problem. Similarly, the SVD equation for a_j (the j^th column of X) is:

$${\mathit{a}}_{\mathit{j}}=\sum _{\mathit{k}=\mathbf{1}}^{\mathit{r}}{\mathit{v}}_{\mathit{jk}}{\mathit{s}}_{\mathit{k}}{\mathit{u}}_{\mathit{k}}$$E4

which is a linear combination of the left singular values {u_k}. The j^th column of V^T, a’_j (see Figure 1), contains the coordinates of the j^th column of X in the coordinate system (basis) of the scaled left singular vectors (the scores), s_ku_k. By using the vector a’_j, the analysis may be captured by r ≤ n variables, which is always fewer than the m elements in the vector a_j, thus SVD reduces the number of variables required. Essentially, there are only r (which we can truncate to eliminate small singular values and further reduce the dimensionality) component vectors (the corresponding right singular vectors) which explain the behavior of X. The application of PCA often use the SVD property:

$${\mathit{X}}^{\left(\mathit{l}\right)}=\sum _{\mathit{k}=\mathbf{1}}^{\mathit{l}}{\mathit{u}}_{\mathit{k}}{\mathit{s}}_{\mathit{k}}{\mathit{v}}_{\mathit{k}}^{\mathit{T}}$$E5

where X^(l) is the closest rank-l matrix to X, i.e., X^(l) minimizes the sum of the squares of the difference between the elements of X and X^(l), diff = ∑_ij|x_ij – x^(l)_ij|².

We can define the covariance matrix as X^TX = Σ_ig_ig_i^T. SVD analysis of X^TX yields V^T, which contains the principal components of {g_i}, i.e. the right singular vectors {v_k} are the same as the principal components of {g_i}. The eigenvalues of X^TX are equivalent to s_k², which are proportional to the variances of the principal components. The matrix XX^T = Σ_ja_ja_j^T is proportional to the covariance matrix of the variables of a_j. The left singular vectors {u_k} are the same as the principal components of {a_j}. The s_k² are proportional to the variances of the principal components. The diagonal values of S (i.e., s_k) are the “singular value spectrum”. The value of a singular value is indicative of its importance in explaining the data. More specifically, the square of each singular value is proportional to the variance explained by each singular vector.

In PCA, X is defined and the SVD is computed. The singular values s_k are then plotted versus k. A reduced dimensionality approximation to X is computed by truncating the singular value series, i.e. by setting s_k = 0 for k > K where K is the chosen threshold value, and using Eq. (3). Note that the squared singular values s²_k are a measure of the variability or “power” in the corresponding signal component specified by the corresponding to singular vector (like the frequency component in Fourier analysis).

2.2 Least squares functional modeling with PCA optimization

We ultimately seek a model that relates the inputs (navigation parameters) to the height error. The least squares method is used to model DEM errors. More details on the least squares method can also be found in Ghilani [17]. Assuming a linear forward model:

$$\mathit{Xh}=e$$E6

The set of navigation data parameters is used to form the h vector. We will assume that the corresponding DEM height errors form the e vector. For each epoch, we form a single h and e vector. Since the model parameters vary with time, vectors are created by stacking the values as a function of time.

Using the matrix form, the X matrix in Eq. (6) represents the mapping between the model parameters and the height. This mapping is what we want to estimate from the DEM error training set. While there are several approaches, the following approach is attractive. Since Eq. (6) must apply for all realizations we create matrices H and E by combination of all the vectors of h and e, i.e., H = [h₁|h₂| … h_n] and similarly for E. We can then write:

$$\mathit{XH}=\mathit{E}$$E7

Shifting the mean and scaling the vector elements affects the performance of the estimate. Assuming we have enough realizations, we can write a least-squares empirical estimate of X, X_e, as:

$${\mathit{X}}_e=\mathit{E}{\mathit{{\rm H}}}^T{\left(\mathit{H}{\mathit{H}}^T\right)}^{-1}$$E8

While we can use X_e directly, SVD analysis of X_e provides us with a powerful tool to understand the relationship between H and E. The SVD analysis of X is:

$${\mathit{X}}_e=\mathit{U\Sigma }{\mathit{V}}^{\mathit{T}}$$E9

The singular value spectrum of Σ tells us the dimensionality (say l) of the model parameter vector (h) are “useful” while the first l column vectors of V provide a basis for the space of “useful” subspace of model parameters. (In the ideal case, a retained column vector of V that contained a single 1 in a particular location w and zero elsewhere would suggest that only the w^th component of model vector needs to be used. In general, the l useful columns of V provide a mapping (rotation and scaling) from the original model parameter space to a restricted model parameter space. As noted previously, the value of l is either the rank or is selected by truncating the singular value spectrum. The first l columns of U form a basis of the output or the range space of X_e, that is, of the height error function. The columns of U and V greater than l can be discarded, which simplifies usage of the truncated model.

Note that while the order of the model parameters or error values within the vectors is not important, the range of cases included does matter. The model may be unable to adequately represent a case not included in the training set. Also, the model can only represent linear functions of the input model parameters. (Recall, that we can produce non-linear responses by including non-linear transforms of model parameters in the model parameter vector.)

Let X_r be the singular value reduced estimate of X_e, i.e. Eq. (9) where the diagonal elements of Σ are set to zero beyond the l^th element. To apply a height correction, we form an h vector of the model parameters from a new data take, and compute an estimate of the height error vector e:

$$X_{f}h=e$$E10

The error e is subtracted from the height map to remove the height error. Note that due to our formulation of the original model and height vectors that this is done “block-wise”.

In this research, the sensitivity analysis is carried out using PCA and the modeling is investigated the least squares function modeling where the system navigation data are considered the inputs to the model and the DEM height error is considered the desired output of the model.

3. Data and methodology

The DEM model errors associated with airborne mapping system datasets are investigated in this research. Figure 2 shows an example of DEM errors case study. It can be seen that the DEM errors map shows large errors in the left side associated with low depression angles and small errors in the right side associated with high depression angles. The across-track DEM errors sections at three different depression angles are investigated (sections 1 to 3) as shown in Figure 2 where the DEM errors are the highest in these locations.

Figure 2.

DEM error map along with the locations of the across-track sections in purple solid lines (1, 2, and 3).

The methodology for the proposed combined PCA analysis and least squares method is shown in Figure 3. The PCA analysis is utilized to identify the significant inputs from multiple navigation data and the least squares method is implemented to estimate the DEM errors models. The root mean squares (RMS) errors are used to quantify the accuracy of the developed DEM errors models.

Figure 3.

Combined PCA analysis and least squares methodology.

4. Results and discussion

Three across-track sections (1 to 3) as shown in Figure 2 have been investigated to test the performance of combined PCA analysis and least squares method. The navigation data (3 positions, 3 velocities, 3 attitudes, 3 accelerations, 3 attitude rates, 3 attitudes accelerations, 1 doppler, and 1 doppler rate) represent the input parameters for DEM errors modeling were investigated using PCA analysis method. Out of 20 inputs, 10 inputs were found significant (2 positions, 2 accelerations, 2 attitudes, 2 attitude rates, 1 doppler and 1 doppler rate) in the three across-track sections.

Then the least squares were employed to model the functional relationship between the 10 significant navigation inputs and the targeted DEM error output. The RMS errors were estimated before and after modeling to test the performance of the developed model to compensate the DEM errors. Figure 4 shows the targeted DEM errors and the modeled DEM errors where Figure 5 shows the differences between the targeted and modeled DEM errors for section (1) case. In section (1) case, the estimated RMS error before modeling is 0.59 m and after modeling is 0.44 m which means that about 25% of DEM errors can be compensated in section (1) case. Figure 6 shows the targeted DEM errors and the modeled DEM errors where Figure 7 shows the differences between the targeted and modeled DEM errors for section (2) case. In section (2) case, the estimated RMS error before modeling is 0.45 m and after modeling is 0.33 m which means that about 27% of DEM errors can be compensated in section (2) case. Figure 8 shows the targeted DEM errors and the modeled DEM errors where Figure 9 shows the differences between the targeted and modeled DEM errors for section (3) case. In section (3) case, the estimated RMS error before modeling is 0.35 m and after modeling is 0.25 m which means that about 28% of DEM errors can be compensated in section (3) case.

Figure 4.

Section (1) DEM errors where blue line represents the targeted DEM error and red lines represents the modeled DEM error.

Figure 5.

Section (1) difference between the targeted DEM error and modeled DEM errors.

Figure 6.

Section (2) DEM errors where blue line represents the targeted DEM error and red lines represents the modeled DEM error.

Figure 7.

Section (2) difference between the targeted DEM error and modeled DEM errors.

Figure 8.

Section (3) DEM errors where blue line represents the targeted DEM error and red lines represents the modeled DEM error.

Figure 9.

Section (3) difference between the targeted DEM error and modeled DEM errors.

Table 1 and Figure 10 summarize the RMS errors and the percentage of DEM error compensation from the three sections where the overall percentage of DEM error compensation is about 27% on average. The results show that the combined PCA analysis and least squares method can be used as a powerful tool to compensate the DEM error produced by the navigation data with about 27% in average for the case study investigated in this research.

Section	RMS error before modeling (m)	RMS error after modeling	Percentage of DEM error compensation (%)
Section (1)	0.59	0.44	25
Section (2)	0.45	0.33	27
Section (3)	0.35	0.25	28
Overall percentage of DEM error compensation			27

Table 1.

RMS errors for all sections before and after DEM error modeling.

Figure 10.

RMS errors for all sections before and after DEM error modeling.

5. Conclusions and recommendation

The sensitivity analysis of DEM errors to the GPS/inertial navigation data was investigated in this research. It was concluded that the sensitivity analysis of the DEM errors can be performed using the PCA to identify the significant GPS/inertial navigation data components that are strongly correlated with DEM errors. Also, it is concluded that the least squares method can be rigorously utilized to establish the functional relationship between the DEM errors and the significant GPS/inertial navigation data components. The combined PCA and least squares method were validated using the DEM model errors associated with airborne mapping system datasets using three across-track sections of the DEM errors datasets. The results show that the combined PCA analysis and least squares method can be used as a powerful tool to compensate for the DEM error produced by the GPS/inertial navigation data with about 27% on average. Therefore, it is recommended to use the combined PCA analysis and least squares method to reduce the DEM errors associated with the DEM model produced by the direct georeferenced airborne mapping system.

Conflict of interest

The author declares no conflict of interest.