# Terrain Referenced Navigation Using a Multilayer Radial Basis Function-Based Extreme Learning Machine

Authors: Jungshin Lee; Changky Sung; Juhyun Oh
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

A high-resolution digital elevation model (DEM) is an important element that determines the performance of terrain referenced navigation (TRN). However, the higher the resolution of the DEM, the bigger the memory size needed for storing it. It is difficult to secure such large memory spaces in small, low-priced unmanned aerial vehicles. In this study, a high-precision terrain regression model to fit the DEM is generated using the extreme learning machine technique based on the multilayer radial basis function. The TRN results using the proposed method are compared with existing studies on various DEM fitting methods. This study verifies that the proposed method obtains improved fitting accuracy and TRN performance over existing DEM fitting methods such as bilinear interpolation, SVM for regression, and bi-spline neural network, without the DEM storage space.


---

## Body

## 1. Introduction

Terrain referenced navigation (TRN) is a navigation technology that estimates the position of an aircraft by comparing the altitude measured by an altimeter with the digital elevation model (DEM). Currently, the global navigation satellite system (GNSS) is commonly used for precise navigation [1] but an alternative is needed for environments in which the GNSS is unavailable due to interference such as hostile jamming. Thus, studies on TRN are actively being carried out as an alternative that will ensure precise navigation performance when GNSS is denied [2–5].The DEM is a map that stores the elevation values of terrain to represent a terrain’s surface in three dimensions, and a high-resolution DEM is essential for improving the performance of TRN. However, large amounts of memory space are needed to store a high-resolution DEM and the loading time of the DEM in a separate storage space interferes with real-time computation. In addition, various interpolation methods are being used to obtain terrain elevation at the desired location using a discontinuous DEM. This method also generates additional computational errors, which can be larger as the resolution of the DEM decreases. To alleviate this problem, studies are continuously being conducted on devising a continuous terrain elevation model through a regression method based on the machine learning technique [6–10]. However, studies on DEM fitting using machine learning techniques up until now were conducted on DEM data in small areas of several hundred square meters or were unable to acquire higher accuracy than interpolation methods. Furthermore, excluding the method proposed by Liu et al. [6], there have been no studies that applied actual TRN using trained regression models. Also, in the case of Liu et al. [6], applying low-resolution DEM learning results in narrow areas would be inappropriate for actual high-performance TRN. In addition to the necessity of very high fitting accuracy for application in TRN, there is also the constraint of guaranteeing real-time operation. Taking into consideration such conditions, we applied the extreme learning machine (ELM) technique to train a DEM. The ELM has a simple structure and a fast learning speed, while it is still able to obtain high levels of accuracy and it is therefore currently being researched in various fields [11–14]. In addition, studies are actively being conducted on advanced algorithms to improve the performance of the conventional ELM or to resolve the complex problems beyond the abilities of the conventional ELM [15–22]. Kan et al. [9] and Hu and Yuan [10] estimated the DEM using a conventional ELM, but it was conducted on a very narrow area, and it was unable to obtain higher accuracy than the Delaunay triangulation-based interpolation method.Even with the advantage of not needing large memory space, unless it has high levels of fitting accuracy that can replace the DEM, it cannot be applied to TRN. Moreover, when fitting a DEM as a complex model to enhance accuracy, it can preclude the real-time operation for navigation. Considering such restrictions, we used the generalized ELM autoencoder [15, 20] and multilayer radial basis function (RBF) [16, 18] to design a model that enhances the fitting accuracy of the ELM and guarantees real-time computing. We compared systems that used the proposed model and bilinear interpolation, which demonstrated the highest TRN performance in past studies, and our results verified that the proposed technique is the most suitable for TRN.In the next section, we introduce existing studies on DEM fitting methods, namely, bilinear interpolation, the bi-spline neural network (B-spline NN), and the support vector machine for regression (SVMR). Then we introduce DEM learning techniques that use the conventional ELM and propose the multilayer RBF-based ELM (ML-RBF-ELM) to overcome the limitations of the conventional ELM. Furthermore, we compare the existing various DEM fitting methods and use of the proposed ML-RBF-ELM regression method and bilinear interpolation. Finally, we verify the design of the proposed technique by comparing it with the TRN performance.


## 2. Conventional DEM Fitting Method

There are different types of DEM database: the shuttle radar topography mission (SRTM) elevation database, digital terrain elevation data (DTED), etc. In this study, we used 10 m resolution DTED level 3 and 30 m resolution DTED level 2. There are various interpolation methods and machine learning regression techniques that can be used to obtain terrain elevation at the desired position using discontinuous DEM data. This section introduces previous studies that are aimed at obtaining continuous terrain elevation information.


### 2.1. Bilinear Interpolation

Interpolation methods frequently used to obtain continuous elevation information include Delaunay triangulation [9] and bilinear and bicubic interpolations. This study introduces bilinear interpolation, which is generally used in TRN systems. Bilinear interpolation is a method for linearly determining the linear distance of the four nearby DEM grid points, Q11x1,y1, Q12x1,y2, Q21x2,y1, and Q22x2,y2, when estimating the elevation value at any point Px,y within the four points as shown in Figure 1.Figure 1
Bilinear interpolation scheme.If the terrain elevation value estimated at the pointPx,y is fx,y, the elevation values stored in DTED, hQ11, hQ12, hQ21, and hQ22, can be used to compute as follows.
(1)fx,y=1x2−x1y2−y1x2−xx−x1hQ11hQ12hQ21hQ22y2−yy−y1.


### 2.2. B-Spline Neural Network

As discussed by Liu et al. [6], the B-spline NN-based regression method was proposed to fit the 30 m resolution DEM. The method divides all the data into multiple bins to supplement the weaknesses of polynomial regression, and it fits the divided data with separate polynomial functions. The two adjacent polynomials’ first and second differentials were made the same for smooth connection without making noncontinuous points in bins where division occurs. Liu et al. [6] designed the system as a single hidden layer where such spline basis functions are used for the layer’s nodes and DEM data were learned as a neural network based on this. In this study, the B-spline NN was designed as 1971 hidden nodes using the spline basis function of order 3. This is compared with other fitting methods in the next section.


### 2.3. Support Vector Machine for Regression

Okwuashi and Ndehedehe [7] used SVMR to fit the DEM in the 225 m by 225 m field. The regression method using SVMR was compared with the artificial neural network and nearest neighbour method, and it showed the highest level of accuracy. SVMR uses the same principles as the SVM for classification with only a few minor differences. For regression, a margin of tolerance was set as the estimate of the SVMR, which would have already been determined from the problem [23]. SVMR should consider more design factors than SVM for classification. However, like SVM for classification, SVMR is learned to minimize the difference between the true value and the predicted value and maximize the margin between decision boundaries. For nonlinear SVMR, kernel functions map input data to higher dimensional feature space, simplifying it as a linear SVMR [23]. In our study, the Gaussian kernel function was used for the precision DTED regression model based on SVMR. In addition, the box constraint was used and set as 0.002178 to prevent overfitting. The optimization issue of the SVMR model is usually solved with the Lagrange dual formulation for computational convenience, and when the Karush-Kuhn-Tucker (KKT) complementarity conditions are satisfied, the optimal solution can be obtained. In this study, such KKT tolerance was set as 0.5×10−4 and it was optimized using the iterative single data algorithm (ISDA), which updates a single multiplier per iteration.


## 2.1. Bilinear Interpolation

Interpolation methods frequently used to obtain continuous elevation information include Delaunay triangulation [9] and bilinear and bicubic interpolations. This study introduces bilinear interpolation, which is generally used in TRN systems. Bilinear interpolation is a method for linearly determining the linear distance of the four nearby DEM grid points, Q11x1,y1, Q12x1,y2, Q21x2,y1, and Q22x2,y2, when estimating the elevation value at any point Px,y within the four points as shown in Figure 1.Figure 1
Bilinear interpolation scheme.If the terrain elevation value estimated at the pointPx,y is fx,y, the elevation values stored in DTED, hQ11, hQ12, hQ21, and hQ22, can be used to compute as follows.
(1)fx,y=1x2−x1y2−y1x2−xx−x1hQ11hQ12hQ21hQ22y2−yy−y1.


## 2.2. B-Spline Neural Network

As discussed by Liu et al. [6], the B-spline NN-based regression method was proposed to fit the 30 m resolution DEM. The method divides all the data into multiple bins to supplement the weaknesses of polynomial regression, and it fits the divided data with separate polynomial functions. The two adjacent polynomials’ first and second differentials were made the same for smooth connection without making noncontinuous points in bins where division occurs. Liu et al. [6] designed the system as a single hidden layer where such spline basis functions are used for the layer’s nodes and DEM data were learned as a neural network based on this. In this study, the B-spline NN was designed as 1971 hidden nodes using the spline basis function of order 3. This is compared with other fitting methods in the next section.


## 2.3. Support Vector Machine for Regression

Okwuashi and Ndehedehe [7] used SVMR to fit the DEM in the 225 m by 225 m field. The regression method using SVMR was compared with the artificial neural network and nearest neighbour method, and it showed the highest level of accuracy. SVMR uses the same principles as the SVM for classification with only a few minor differences. For regression, a margin of tolerance was set as the estimate of the SVMR, which would have already been determined from the problem [23]. SVMR should consider more design factors than SVM for classification. However, like SVM for classification, SVMR is learned to minimize the difference between the true value and the predicted value and maximize the margin between decision boundaries. For nonlinear SVMR, kernel functions map input data to higher dimensional feature space, simplifying it as a linear SVMR [23]. In our study, the Gaussian kernel function was used for the precision DTED regression model based on SVMR. In addition, the box constraint was used and set as 0.002178 to prevent overfitting. The optimization issue of the SVMR model is usually solved with the Lagrange dual formulation for computational convenience, and when the Karush-Kuhn-Tucker (KKT) complementarity conditions are satisfied, the optimal solution can be obtained. In this study, such KKT tolerance was set as 0.5×10−4 and it was optimized using the iterative single data algorithm (ISDA), which updates a single multiplier per iteration.


## 3. DEM Fitting Method by ML-RBF-ELM

### 3.1. Conventional ELM

Since first being developed by Huang et al., the ELM has demonstrated good generalization ability and much faster learning than the backpropagation method of stochastic gradient descent and is therefore being studied in various fields. The ELM is a single hidden layer feedforward neural network. The input weights and hidden unit biases are chosen randomly and do not need to be adjusted [24]. As the output weights are analytically computed by the least squares method, this method provides a good generalization performance at an extremely fast learning speed. For M distinct training sample ℵ=xm,dmxm∈RN,dm∈RL,m=1,⋯,M, the output of ELM with hidden layer nodes is as in [8, 9]:
(2)ϕk=exp−βkxim−aki2,k=1,⋯,K,i=1,⋯,N.Here,ϕk is the activation function of the kth hidden node, which is modelled as a RBF. xm=xim,i=1,⋯,N is the mth input data and N is set to 2 in this study. x1m and x2m are the normalized longitude and latitude, respectively. dm is the mth target (or true) data and means the normalized terrain height of the DEM in this study. αki and βk are the center and width, respectively, of the kth RBF. By Huang et al. [24], the center and width are assigned randomly. ∙ means the Euclidean distance. K is the number of hidden nodes and is set to 271 including a bias node in this study. The number of nodes was optimized to have the lowest training and test errors. The ELM can estimate any target function with zero error using equation (2). The relation is written as
(3)HC=T,H=hx1⋮hxM=ϕ1α1,β1,x1⋯ϕKαK,βK,x11⋮⋱⋮⋮ϕ1α1,β1,xM⋯ϕKαK,βK,xM1.Here,HM×K is the hidden layer output matrix and hxm is the hidden layer set of the mth sample data. The hidden layer set consists of K activation functions. TM×L=d1,⋯,dm,⋯,dMT is the target matrix of the ELM, and L is set to 1 because the normalized terrain elevation data is used as the true value of the training data. CK×L is the weight matrix of the output layer and the smallest norm least squares solution of equation (3) as follows:
(4)C=H†T.Here,H† is the Moore-Penrose generalized inverse, and it is generally calculated as a singular value composition. In order to prevent the overfitting and improve the robustness, a ridge parameter is used [15, 16]. The optimization problem is described as
(5)min12CTC+12λ∑m=1Mξm2.Here,λ is a regularization parameter and is set to 0.2. ξm is the residual between the target value dm and estimated output value ym of the mth sample. According to KKT conditions, the solution of equation (5) is calculated as follows [9]:
(6)C=H†T=HTH+λI−1HTT,M≥K,HTHHT+λI−1T,M<K.Here,I is a L-by-L identity matrix.


## 3.1. Conventional ELM

Since first being developed by Huang et al., the ELM has demonstrated good generalization ability and much faster learning than the backpropagation method of stochastic gradient descent and is therefore being studied in various fields. The ELM is a single hidden layer feedforward neural network. The input weights and hidden unit biases are chosen randomly and do not need to be adjusted [24]. As the output weights are analytically computed by the least squares method, this method provides a good generalization performance at an extremely fast learning speed. For M distinct training sample ℵ=xm,dmxm∈RN,dm∈RL,m=1,⋯,M, the output of ELM with hidden layer nodes is as in [8, 9]:
(2)ϕk=exp−βkxim−aki2,k=1,⋯,K,i=1,⋯,N.Here,ϕk is the activation function of the kth hidden node, which is modelled as a RBF. xm=xim,i=1,⋯,N is the mth input data and N is set to 2 in this study. x1m and x2m are the normalized longitude and latitude, respectively. dm is the mth target (or true) data and means the normalized terrain height of the DEM in this study. αki and βk are the center and width, respectively, of the kth RBF. By Huang et al. [24], the center and width are assigned randomly. ∙ means the Euclidean distance. K is the number of hidden nodes and is set to 271 including a bias node in this study. The number of nodes was optimized to have the lowest training and test errors. The ELM can estimate any target function with zero error using equation (2). The relation is written as
(3)HC=T,H=hx1⋮hxM=ϕ1α1,β1,x1⋯ϕKαK,βK,x11⋮⋱⋮⋮ϕ1α1,β1,xM⋯ϕKαK,βK,xM1.Here,HM×K is the hidden layer output matrix and hxm is the hidden layer set of the mth sample data. The hidden layer set consists of K activation functions. TM×L=d1,⋯,dm,⋯,dMT is the target matrix of the ELM, and L is set to 1 because the normalized terrain elevation data is used as the true value of the training data. CK×L is the weight matrix of the output layer and the smallest norm least squares solution of equation (3) as follows:
(4)C=H†T.Here,H† is the Moore-Penrose generalized inverse, and it is generally calculated as a singular value composition. In order to prevent the overfitting and improve the robustness, a ridge parameter is used [15, 16]. The optimization problem is described as
(5)min12CTC+12λ∑m=1Mξm2.Here,λ is a regularization parameter and is set to 0.2. ξm is the residual between the target value dm and estimated output value ym of the mth sample. According to KKT conditions, the solution of equation (5) is calculated as follows [9]:
(6)C=H†T=HTH+λI−1HTT,M≥K,HTHHT+λI−1T,M<K.Here,I is a L-by-L identity matrix.


## 4. ML-RBF-ELM

The conventional ELM is easy to configure and the learning speed is fast, but it cannot be applied to fields that process big data with very large nonlinearity or require high levels of accuracy. To solve these problems, many studies are being actively conducted on various advanced ELM algorithms. In the constraint in which real-time computation must be guaranteed, three hidden layers were used as shown in Figure2 to maximize the training and test accuracies. The first hidden layer was designed as an ELM-based autoencoder (ELM-AE), and the remaining hidden layers were designed as conventional ELMs. As the main difference between the ELM-AE and conventional ELM, the ELM is a supervised NN and the outputs are the target values but the ELM-AE is an unsupervised NN and the outputs are equal to the inputs [15].Figure 2
Proposed ML-RBF-ELM architecture.The output weight matrix of the 1st hidden layer outputC1 is calculated as follows [15]:
(7)C1=HX1†=H1TH1+λI−1H1TX,M≥K1,H1TH1H1T+λI−1X,M<K1.Here,XM×2=x1,⋯,xm,⋯,xMTand xm=xim,i=1,⋯,N. I is a N-by-N identity matrix and K1 is the number of nodes in the 1st hidden layer. In this study, the ELM-AE was used to make the inputs more important in a rough terrain area than in a flat area. The 2nd and 3rd hidden layers shown in Figure 2 are similar to the conventional ELM. These layers were designed to represent the terrain height with refined accuracy. Unlike the conventional ELM, in the proposed ML-RBF-ELM method, the weights between the previous hidden layer and the current hidden layer were determined using K-means algorithms to enhance the training accuracy. The 1st and 2nd intermediate vectors are virtual vectors that represent the process in which the output weight matrices of the 1st and 2nd hidden layers explained in Algorithm 1 are learned. K-means clustering is used to partition M sample data into K clusters based on features in which each sample data belongs to the cluster with the nearest mean of a Gaussian function [25]. In this study, the mean and width values of the RBF nodes in the hidden layers are determined by the K-mean clustering based on the expectation-maximization algorithm [25]. The process of the proposed ML-RBF-ELM is given in Algorithm 1. The numbers of nodes in the three hidden layers, K1, K2, and K3, are set to 384, 2301, and 193, respectively.Algorithm 1: Proposed ML-RBF-ELM algorithm.
(1) Prepare a training setℵ=xm,dm∣xm∈R2,dm∈R1,m=1,⋯,M where xm=x1m,x2m is the input data with the normalized longitude and latitude and dm is the normalized terrain height in the DEM
(2) Set the activation function of the 1st hidden layer,Φ1=ϕ11,⋯,ϕ1K1, by equation (2) and assign the center α1k1i and width β1k1 of the RBF by K-means clustering
(3) Calculate the 1st hidden layer output matrixH1 by equation (3)
(4) Calculate the output weight matrix of the 1st hidden layer,C1K1+1×2=H1TH1+λI−1H1TX where XM×2=x1,⋯,xMT
(5) Calculate the new input data of the 2nd hidden layer:X′M×2=H1C1
(6) Set the activation function of the 2nd hidden layer,Φ2=ϕ21,⋯,ϕ2K2, by equation (2) and assign the center α2k2i and width β2k2 of the RBF by K-means clustering
(7) Calculate the 2nd hidden layer output matrixH2 by equation (3)
(8) Calculate the output weight matrix of the 2nd hidden layer,C2K2+1×1=H2TH2+λI−1H2TT where TM×1=d1,⋯,dMT
(9) Calculate the new input data of the 3rd hidden layer:y′M×1=H2C2
(10) Set the activation function of the 3rd hidden layer,Φ3=ϕ31,⋯,ϕ3K3, by equation (2) and assign the center α3k3i and width β3k3 of the RBF by K-means clustering
(11) Calculate the 3rd hidden layer output matrixH3 by equation (3)
(12) Calculate the output weight matrix of the 3rd hidden layer,C3K3+1×1=H3TH3+λI−1H3TT
(13) Compute and verify the regression results,yM×1=H3C3


### 4.1. Comparison of Various DEM Fitting Methods

Figure3 illustrates the DEM fitting area with a total area of about 689.70 km2 for DTED level 3 and 4203.86 km2 for DTED level 2. As mentioned above, we used 10 m resolution DTED level 3 and 30 m resolution DTED level 2 for the training and test data. The total obtainable number of sample data were 5391738 sets (161 MB) for level 3 and 5443200 sets (163 MB) for level 2. At a rate of 8 : 2, the sample data was used by dividing the training data (4313390 sets, 129 MB for level 3, and 4354560 sets, 130.4 MB for level 2) and test data (1078348 sets, 32 MB for level 3, and 1088640 sets, 32.6 MB for level 2). The sample data are acquired from the grid points such as Q11, Q12, Q21, and Q22 in Figure 1. The learning was performed with an Intel® Core™ i7-6700K, CPU @4.00 GHz, RAM 16.0 GB computing environment. All simulations for the bilinear interpolation, B-spline NN, SVMR, conventional ELM, and the proposed ML-RBF-ELM were carried out in a MATLAB R2018a environment.Figure 3
DEM fitting areas and flight trajectories. (a) DEM level 3 fitting area. (b) DEM level 2 fitting area.
(a)
(b)Table1 compares the fitting accuracy and operation time according to the bilinear interpolation, B-spline NN, SVMR, conventional ELM, and the proposed ML-RBF-ELM for fitting DTED level 3. The fitting results of level 2 are also similar to those of Table 1. Each method was designed to maximize fitting accuracy. For the ML-RBF-ELM design, which was designed with three hidden layers, as shown in Table 1, the results confirmed that compared to the B-spline NN and conventional ELM designed as a single hidden layer, more nodes could be used in each hidden layer. Furthermore, we also confirmed that the ML-RBF-ELM, which was designed with three hidden layers, learned faster than the B-spline NN, which was designed as a single hidden layer. Likewise, it was evident that the ELM performed better in a shorter amount of time than the NN. The test time in Table 1 indicates the time needed for fitting one new sample to the terrain height. For the index to evaluate the fitting accuracy, the following mean value of the absolute percentage error was used. The value E¯ is calculated as follows:
(8)E¯=∑m=1Mym−dmdm×100×1M.Table 1
Results of various DEM fitting methods.
Bilinear interpolation
B-spline NN
SVMR
Conventional ELM
ML-RBF-ELM
Number of hidden nodes
—
1971
—
271
3842301193
Training time
—
249.3 sec
3135.4 sec
8.10 sec
71.72 sec
Test time
—
99.8μsec
115.3 sec
13.0μsec
44.0μsec
Training error
Approx. 1.02%
3.81%
8.76%
24.03%
0.98%
Test error
Approx. 1.01%
7.50%
8.64%
23.48%
1.10%Here,ym and dm are the estimated and target output value of the mth sample data, respectively. M is the total number of samples.The training and test errors of the bilinear interpolation in Table1 are approximate values. These values are about 25% of the absolute percentage errors when the terrain elevation value estimated at point Q21 was acquired from points Q11, Q31, Q20, and Q22 in Figure 1. The conventional ELM had the smallest training and test time but the lowest fitting accuracy so it is not suitable for TRN. The SVMR did not satisfy the test time restrictions for real-time operation, so it is unsuitable. However, the proposed ML-RBF-ELM had the best fitting accuracy and its short test time allowed for real-time operation, so we judged that it would be suitable in TRN. Figures 4 and 5 are illustrations of the terrain height fitting results for the test data. We found that previous research results did not reach the accuracy of bilinear interpolation, as shown in Figure 4. Despite the advantages of not requiring a large memory space to store the DEM and the time needed for loading, if fitting accuracy results cannot be obtained at the same levels of bilinear interpolation, then the performance of TRN cannot be guaranteed. On the other hand, as evident in Figure 5, the proposed ML-RBF-ELM demonstrates a fitting accuracy similar to the bilinear interpolation of Figure 4. To store the weights of the proposed ML-RBF-ELM, we need a 105 KB size memory. This is only about 1/1533 of the 161 MB memory size for DTED level 3 and about 1/1552 of the 163 MB memory size for DTED level 2.Figure 4
Results of the DEM fitting methods using DTED level 3. (a) Bilinear interpolation. (b) B-spline neural network. (c) SVMR. (d) Conventional ELM.
(a)
(b)
(c)
(d)Figure 5
Fitting result of the proposed ML-RBF-ELM. (a) DTED level 3 result. (b) DTED level 2 result.
(a)
(b)


## 4.1. Comparison of Various DEM Fitting Methods

Figure3 illustrates the DEM fitting area with a total area of about 689.70 km2 for DTED level 3 and 4203.86 km2 for DTED level 2. As mentioned above, we used 10 m resolution DTED level 3 and 30 m resolution DTED level 2 for the training and test data. The total obtainable number of sample data were 5391738 sets (161 MB) for level 3 and 5443200 sets (163 MB) for level 2. At a rate of 8 : 2, the sample data was used by dividing the training data (4313390 sets, 129 MB for level 3, and 4354560 sets, 130.4 MB for level 2) and test data (1078348 sets, 32 MB for level 3, and 1088640 sets, 32.6 MB for level 2). The sample data are acquired from the grid points such as Q11, Q12, Q21, and Q22 in Figure 1. The learning was performed with an Intel® Core™ i7-6700K, CPU @4.00 GHz, RAM 16.0 GB computing environment. All simulations for the bilinear interpolation, B-spline NN, SVMR, conventional ELM, and the proposed ML-RBF-ELM were carried out in a MATLAB R2018a environment.Figure 3
DEM fitting areas and flight trajectories. (a) DEM level 3 fitting area. (b) DEM level 2 fitting area.
(a)
(b)Table1 compares the fitting accuracy and operation time according to the bilinear interpolation, B-spline NN, SVMR, conventional ELM, and the proposed ML-RBF-ELM for fitting DTED level 3. The fitting results of level 2 are also similar to those of Table 1. Each method was designed to maximize fitting accuracy. For the ML-RBF-ELM design, which was designed with three hidden layers, as shown in Table 1, the results confirmed that compared to the B-spline NN and conventional ELM designed as a single hidden layer, more nodes could be used in each hidden layer. Furthermore, we also confirmed that the ML-RBF-ELM, which was designed with three hidden layers, learned faster than the B-spline NN, which was designed as a single hidden layer. Likewise, it was evident that the ELM performed better in a shorter amount of time than the NN. The test time in Table 1 indicates the time needed for fitting one new sample to the terrain height. For the index to evaluate the fitting accuracy, the following mean value of the absolute percentage error was used. The value E¯ is calculated as follows:
(8)E¯=∑m=1Mym−dmdm×100×1M.Table 1
Results of various DEM fitting methods.
Bilinear interpolation
B-spline NN
SVMR
Conventional ELM
ML-RBF-ELM
Number of hidden nodes
—
1971
—
271
3842301193
Training time
—
249.3 sec
3135.4 sec
8.10 sec
71.72 sec
Test time
—
99.8μsec
115.3 sec
13.0μsec
44.0μsec
Training error
Approx. 1.02%
3.81%
8.76%
24.03%
0.98%
Test error
Approx. 1.01%
7.50%
8.64%
23.48%
1.10%Here,ym and dm are the estimated and target output value of the mth sample data, respectively. M is the total number of samples.The training and test errors of the bilinear interpolation in Table1 are approximate values. These values are about 25% of the absolute percentage errors when the terrain elevation value estimated at point Q21 was acquired from points Q11, Q31, Q20, and Q22 in Figure 1. The conventional ELM had the smallest training and test time but the lowest fitting accuracy so it is not suitable for TRN. The SVMR did not satisfy the test time restrictions for real-time operation, so it is unsuitable. However, the proposed ML-RBF-ELM had the best fitting accuracy and its short test time allowed for real-time operation, so we judged that it would be suitable in TRN. Figures 4 and 5 are illustrations of the terrain height fitting results for the test data. We found that previous research results did not reach the accuracy of bilinear interpolation, as shown in Figure 4. Despite the advantages of not requiring a large memory space to store the DEM and the time needed for loading, if fitting accuracy results cannot be obtained at the same levels of bilinear interpolation, then the performance of TRN cannot be guaranteed. On the other hand, as evident in Figure 5, the proposed ML-RBF-ELM demonstrates a fitting accuracy similar to the bilinear interpolation of Figure 4. To store the weights of the proposed ML-RBF-ELM, we need a 105 KB size memory. This is only about 1/1533 of the 161 MB memory size for DTED level 3 and about 1/1552 of the 163 MB memory size for DTED level 2.Figure 4
Results of the DEM fitting methods using DTED level 3. (a) Bilinear interpolation. (b) B-spline neural network. (c) SVMR. (d) Conventional ELM.
(a)
(b)
(c)
(d)Figure 5
Fitting result of the proposed ML-RBF-ELM. (a) DTED level 3 result. (b) DTED level 2 result.
(a)
(b)


## 5. Performance of TRN Using ML-RBF-ELM

### 5.1. BKF-Based TRN and INS/TRN Integration

We verified the ML-RBF-ELM by conducting TRN with the ML-RBF-ELM and bilinear interpolation, which is superior to the fitting methods explained in the previous section. TRN simulations were performed using the interferometric radar altimeter (IRA) measurements and a bank of Kalman filters (BKF), which were used in Heli/SITAN [26]. Figure 6 presents a schematic diagram of our system with an inertial navigation system (INS), IRA, barometer, DEM, TRN S/W, and INS/TRN S/W. The TRN progresses in two steps, the time propagation and the measurement update. If the IRA measurements do not satisfy validity check conditions; the TRN performs only the time propagation. To acquire all the navigation information, including the position, velocity, and attitude, the INS/TRN-integrated navigation is needed. Therefore, the loosely coupled INS/TRN-integrated navigation was designed with an extended Kalman filter (EKF) and uses the position estimated by the BKF-based TRN for measurement. Details of the INS/TRN-integrated navigation are given in section 4.3 of reference [27]. If the estimates of TRN do not satisfy terrain validity check conditions, the INS/TRN performs only the time propagation. Table 2 shows the simulation conditions and BKF design parameters. On June 21, 2018, and January 21, 2019, two captive flight tests were carried out on a CESSNA 172 Skyhawk. The INS/GNSS-integrated navigation data and IRA outputs were used as the true trajectory and measurements in this simulation, respectively. We conducted Monte Carlo simulations 100 times using bilinear interpolation and the proposed ML-RBF-ELM. The proposed ML-RBF-ELM was implemented by using a S/W without the DEM, as shown in Figure 6.Figure 6
Schematic diagram of the TRN and INS/TRN system.Table 2
Simulation conditions and BKF design parameters.
Parameter
Value
Initial covariance
152 m2
Accelerometer bias
100μg
Gyro bias
0.01 deg/hr
Gyro white noise
0.005 deg/hr
Barometer bias
14 m
Barometer scale factor
0.2% of height
Barometer white noise
5 m
Process noise covariance
42 m2
Measurement noise covariance
302 m2
Update frequency
50 HzThe BKF is a linear filter that uses the multiple-model adaptive estimation (MMSE) technique to acquire navigation information sequentially, even in the case of a large initial position error. A state variable of one-dimensional Kalman filters comprising the BKF is the vertical bias error, which gradually changes. The BKF forms grids that are centered on the INS/TRN position estimate and assigns a Kalman filter to each grid. In this study,100×100 filters in the grid region corresponding to three times the initial position error were allocated. The time propagation for the jth (1<j≤100×100) filter state variable xjk and its covariance estimate pjk assigned to each grid is described as
(9)xjk=xjk−1+uk,pjk=pjk−1+Qk∆t.Here,Qk is the process noise covariance and ∆t is the sampling time. uk is the upward velocity of the aircraft in the current step. In this study, the IRA was used to measure the angle perpendicular to the direction of flight and the range from the aircraft to the nearest terrain. It then converted these measurements to three-dimension position information. As shown in Figure 7, the relative position vector δxIRA of the nearest point from the aircraft is calculated as follows:
(10)δxIRA=δλresδϕreshres=ρcosξsinαsinβ+ρsinξcosβρcosξsinαcosβ−ρsinξsinβρcosξcosα.Figure 7
Relative position to the nearest point.Here,δλres, δϕres, and hres are the relative distances in the directions of longitude, latitude, and height, respectively. ρ and ξ are the range and look angle output, respectively, of the IRA. The virtual pitch angle α and azimuth angle β of the zero Doppler line are determined by the velocity of the aircraft [27].Thejth filter gain kjk and the measured values zjk using the IRA measurement δxIRA are as follows:
(11)kjk=pjk∙pjk+Rk−1,zjk=hjkλ^k+δλresRew−1,ϕ^k+δϕresRns−1+hres−hbaro.Here,λ^k and ϕ^k are the longitude and latitude, respectively, estimated by INS/TRN and Rk is the measurement noise covariance. λ^k+δλresRew−1 and ϕ^k+δϕresRns−1 are the longitude and latitude, respectively, at the nearest point T as shown in Figure 7. Rew and Rns represent the radius of the curvature of the Earth ellipsoid in the north-south and east-west, respectively. hjk is the terrain height of the jth filter that is centered on the nearest point T. The terrain height is calculated using the bilinear interpolation from the stored DEM or the regression model by the proposed ML-RBF-ELM. hbaro is the aircraft altitude measured by the barometer. The weighted residual squares (WRS) of the jth filter at the kth step WRSjk are calculated as follows [26]:
(12)WRSjk=zjk−xjk2∙pjk+Rk−1.The WRS indicates how well each filter in the grids matches the designed model. The smoothed WRS (SWRS) is expressed in the same way as the moving average of the WRS as follows [26]:
(13)SWRSjk=αWRSjk+1−αSWRSjk−1.Here,SWRSj0=1 and α is the smoothing constant (0.0<α≤1.0) and is determined by calculating how fast the INS estimate gets out of the 3×3 filter exclusion region by the INS navigation error. In this study, α was set as 0.018. The estimated covariance NSWRSi of each SWRS value in the 3×3 filter exclusion region centered on the grid with the minimum SWRS value SWRSmin among total grids is as follows:
(14)NWRSi=SWRSiSWRSmin−1,1≤i≤9.Using the estimated covariance, the longitude and latitudexλk and xϕk estimated by the BKF are calculated as follows:
(15)xλk=xλ0+ve∆t+peRew−1,xϕk=xϕ0+vn∆t+pnRns−1,pi=exp−0.5NWRSi∙∑i=19exp−0.5NWRSi−1,pe=∑i=19peipi,pn=∑i=19pnipi.Here,xλ0 and xϕ0 are the initial longitude and latitude, respectively. pei and pni represent the east and north axis position estimates of the ith filter in the exclusion region. The terrain validity check conditions in Figure 6 are as follows:
(16)SWRSmin∗−SWRSminSWRSmin−1>thr,NMS>thu.Here,SWRSmin∗ is the minimum SWRS value outside the beta area and NMS represents how many times the filter retains the minimum SWRS value. thr and thu are their respective thresholds. Since the smallest value in the grid is SWRSmin, the terrain validity check conditions are always positive and the larger the difference, the rougher the terrain. The more often the above equation is satisfied, the more likely it is that the terrain is unique. The thresholds were each set as 5×10−4 and 2 for the rough area and 5×10−5 and 29 for the flat area through numerical simulations.


## 5.1. BKF-Based TRN and INS/TRN Integration

We verified the ML-RBF-ELM by conducting TRN with the ML-RBF-ELM and bilinear interpolation, which is superior to the fitting methods explained in the previous section. TRN simulations were performed using the interferometric radar altimeter (IRA) measurements and a bank of Kalman filters (BKF), which were used in Heli/SITAN [26]. Figure 6 presents a schematic diagram of our system with an inertial navigation system (INS), IRA, barometer, DEM, TRN S/W, and INS/TRN S/W. The TRN progresses in two steps, the time propagation and the measurement update. If the IRA measurements do not satisfy validity check conditions; the TRN performs only the time propagation. To acquire all the navigation information, including the position, velocity, and attitude, the INS/TRN-integrated navigation is needed. Therefore, the loosely coupled INS/TRN-integrated navigation was designed with an extended Kalman filter (EKF) and uses the position estimated by the BKF-based TRN for measurement. Details of the INS/TRN-integrated navigation are given in section 4.3 of reference [27]. If the estimates of TRN do not satisfy terrain validity check conditions, the INS/TRN performs only the time propagation. Table 2 shows the simulation conditions and BKF design parameters. On June 21, 2018, and January 21, 2019, two captive flight tests were carried out on a CESSNA 172 Skyhawk. The INS/GNSS-integrated navigation data and IRA outputs were used as the true trajectory and measurements in this simulation, respectively. We conducted Monte Carlo simulations 100 times using bilinear interpolation and the proposed ML-RBF-ELM. The proposed ML-RBF-ELM was implemented by using a S/W without the DEM, as shown in Figure 6.Figure 6
Schematic diagram of the TRN and INS/TRN system.Table 2
Simulation conditions and BKF design parameters.
Parameter
Value
Initial covariance
152 m2
Accelerometer bias
100μg
Gyro bias
0.01 deg/hr
Gyro white noise
0.005 deg/hr
Barometer bias
14 m
Barometer scale factor
0.2% of height
Barometer white noise
5 m
Process noise covariance
42 m2
Measurement noise covariance
302 m2
Update frequency
50 HzThe BKF is a linear filter that uses the multiple-model adaptive estimation (MMSE) technique to acquire navigation information sequentially, even in the case of a large initial position error. A state variable of one-dimensional Kalman filters comprising the BKF is the vertical bias error, which gradually changes. The BKF forms grids that are centered on the INS/TRN position estimate and assigns a Kalman filter to each grid. In this study,100×100 filters in the grid region corresponding to three times the initial position error were allocated. The time propagation for the jth (1<j≤100×100) filter state variable xjk and its covariance estimate pjk assigned to each grid is described as
(9)xjk=xjk−1+uk,pjk=pjk−1+Qk∆t.Here,Qk is the process noise covariance and ∆t is the sampling time. uk is the upward velocity of the aircraft in the current step. In this study, the IRA was used to measure the angle perpendicular to the direction of flight and the range from the aircraft to the nearest terrain. It then converted these measurements to three-dimension position information. As shown in Figure 7, the relative position vector δxIRA of the nearest point from the aircraft is calculated as follows:
(10)δxIRA=δλresδϕreshres=ρcosξsinαsinβ+ρsinξcosβρcosξsinαcosβ−ρsinξsinβρcosξcosα.Figure 7
Relative position to the nearest point.Here,δλres, δϕres, and hres are the relative distances in the directions of longitude, latitude, and height, respectively. ρ and ξ are the range and look angle output, respectively, of the IRA. The virtual pitch angle α and azimuth angle β of the zero Doppler line are determined by the velocity of the aircraft [27].Thejth filter gain kjk and the measured values zjk using the IRA measurement δxIRA are as follows:
(11)kjk=pjk∙pjk+Rk−1,zjk=hjkλ^k+δλresRew−1,ϕ^k+δϕresRns−1+hres−hbaro.Here,λ^k and ϕ^k are the longitude and latitude, respectively, estimated by INS/TRN and Rk is the measurement noise covariance. λ^k+δλresRew−1 and ϕ^k+δϕresRns−1 are the longitude and latitude, respectively, at the nearest point T as shown in Figure 7. Rew and Rns represent the radius of the curvature of the Earth ellipsoid in the north-south and east-west, respectively. hjk is the terrain height of the jth filter that is centered on the nearest point T. The terrain height is calculated using the bilinear interpolation from the stored DEM or the regression model by the proposed ML-RBF-ELM. hbaro is the aircraft altitude measured by the barometer. The weighted residual squares (WRS) of the jth filter at the kth step WRSjk are calculated as follows [26]:
(12)WRSjk=zjk−xjk2∙pjk+Rk−1.The WRS indicates how well each filter in the grids matches the designed model. The smoothed WRS (SWRS) is expressed in the same way as the moving average of the WRS as follows [26]:
(13)SWRSjk=αWRSjk+1−αSWRSjk−1.Here,SWRSj0=1 and α is the smoothing constant (0.0<α≤1.0) and is determined by calculating how fast the INS estimate gets out of the 3×3 filter exclusion region by the INS navigation error. In this study, α was set as 0.018. The estimated covariance NSWRSi of each SWRS value in the 3×3 filter exclusion region centered on the grid with the minimum SWRS value SWRSmin among total grids is as follows:
(14)NWRSi=SWRSiSWRSmin−1,1≤i≤9.Using the estimated covariance, the longitude and latitudexλk and xϕk estimated by the BKF are calculated as follows:
(15)xλk=xλ0+ve∆t+peRew−1,xϕk=xϕ0+vn∆t+pnRns−1,pi=exp−0.5NWRSi∙∑i=19exp−0.5NWRSi−1,pe=∑i=19peipi,pn=∑i=19pnipi.Here,xλ0 and xϕ0 are the initial longitude and latitude, respectively. pei and pni represent the east and north axis position estimates of the ith filter in the exclusion region. The terrain validity check conditions in Figure 6 are as follows:
(16)SWRSmin∗−SWRSminSWRSmin−1>thr,NMS>thu.Here,SWRSmin∗ is the minimum SWRS value outside the beta area and NMS represents how many times the filter retains the minimum SWRS value. thr and thu are their respective thresholds. Since the smallest value in the grid is SWRSmin, the terrain validity check conditions are always positive and the larger the difference, the rougher the terrain. The more often the above equation is satisfied, the more likely it is that the terrain is unique. The thresholds were each set as 5×10−4 and 2 for the rough area and 5×10−5 and 29 for the flat area through numerical simulations.


## 6. Simulation Results

Figure8 compares the case of using the bilinear interpolation in the same simulation conditions (dotted lines) and the case of using the regression model fitted according to the proposed ML-RBF-ELM (solid lines). The position root mean square (RMS) errors in the north-south and east-west axes were smaller for the BKF-based TRN using the ML-RBF-ELM. Figure 8(a) shows such characteristics more clearly between 400 seconds and 500 seconds. Figure 8(b) also shows improved performance from 500 seconds onward. The 100 Monte Carlo simulation results of INS/TRN-integrated navigation, which uses the latitude and longitude information estimated by TRN as measurements, are shown in Figure 9. The performance difference between the bilinear interpolation and the ML-RBF-ELM in INS/TRN is smaller compared to TRN. However, it is evident that compared to using bilinear interpolation in Figures 9(a) and 9(c), using ML-RBF-ELM carries out TRN and INS/TRN more stably as shown in Figures 9(b) and 9(d). If the terrain height changes smoothly, the error of bilinear interpolation is small. However, the bilinear interpolation error may be bigger in terrains with higher nonlinearity. Meanwhile, the proposed ML-RBF-ELM is less affected by such terrain roughness or uniqueness and so the change in the overall RMS error may be smaller than bilinear interpolation. Table 3 shows the values of converting the horizontal position RMS errors of TRN and INS/TRN to circular error probability (CEP). Comparing TRN with the bilinear interpolation, TRN with the ML-RBF-ELM improved the result by about 2.335 m CEP in DTED level 3 and about 1.402 m CEP in DTED level 2. In other words, using ML-RBF-ELM, it is possible to fit a highly accurate regression model without the DEM. Moreover, TRN and INS/TRN using the proposed method achieved stable performance.Figure 8
BKF-based TRN RMS errors. (a) DTED level 3 result. (b) DTED level 2 result.
(a)
(b)Figure 9
INS/TRN RMS errors. (a) Bilinear interpolation with DTED level 3. (b) ML-RBF-ELM with DTED level 3. (c) Bilinear interpolation with DTED level 2. (b) ML-RBF-ELM with DTED level 2.
(a)
(b)
(c)
(d)Table 3
Comparison of TRN and INS/TRN RMS errors.
DTED level
RMS errors
Bilinear interpolation
ML-RBF-ELM
3
TRN
7.465 m CEP
5.130 m CEP
3
INS/TRN
4.126 m CEP
3.793 m CEP
2
TRN
12.540 m CEP
11.138 m CEP
2
INS/TRN
9.891 m CEP
6.079 m CEP


## 7. Conclusions

This study introduced DEM fitting methods, a key technology of TRN, that can be used as an alternative in environments where GNSS is unavailable. The DEM is a database that stores the terrain’s elevation data at a constant resolution. Therefore, large amounts of memory are needed to use the high-resolution DEM. Furthermore, the various interpolation methods for obtaining continuous elevation information from the discontinuous DEM includes fitting errors. Accordingly, there are currently many ongoing studies on DEM learning through various machine learning techniques. However, these efforts have only produced results on DEM learning in low-resolution or narrow fields and they have yet to achieve stable performance for TRN using these fitting methods. In this study, we used the ML-RBF-ELM to fit a 10 m resolution DEM in an area of680.70 km2 and a 30 m resolution DEM in an area of 4203.86 km2. The proposed regression model by the ML-RBF-ELM operates more quickly than the previous machine learning regression methods for the DEM and demonstrates similar fitting results with bilinear interpolation. Furthermore, we used BKF-based TRN, which showed improved TRN performance results over using the bilinear interpolation. Thus, our study verified that the proposed ML-RBF-ELM can make real-time operation of navigation and is suitable to TRN. It is clear that the proposed technology can be used by low-priced small unmanned aerial vehicles to execute TRN in environments where GNSS is not possible. Moreover, the proposed ML-RBF-ELM can be used in the field of simultaneous localization and map building (SLAM) because of the small training time.



---

*Source: 9142694-2019-07-11.xml*