Validation of an AI-Assisted Terrain Aided Navigation Algorithm Using Real-World Flight Test Instrumentation Data

1. Introduction

This study presents an AI-enhanced Terrain Aided Navigation (TAN) designed for Hürjet, an advanced jet trainer aircraft seen in Figure 1, with the aim of improving navigation accuracy and operational efficiency. The system utilizes an advanced terrain server combined with a high-resolution digital elevation map (DEM) called digital terrain elevation data (DTED) to produce accurate height-on-terrain (HoT) and height-above-terrain (HaT) data, which are crucial for precise low-altitude aircraft navigation.

The TAN is fundamentally based on a sophisticated technique that includes intricate filter setups. We specifically implement and evaluate banks of adaptive filters, which dynamically react to varying terrain roughness and other environmental factors. This adaptive filtering method enhances the system’s tracking accuracy and considerably improves its resilience under diverse flying situations. Additionally, the TAN incorporates an AI-driven supervisory module that can automatically switch among various software designs, providing a warrant that the system selects the most efficient mode’s output at all times. A simplified diagram of the algorithm is given in Figure 2 where latitude $\delta L$, longitude $\delta \lambda$ , and altitude $\delta h$ errors in geodetic coordinate system. Also, $h_{radalt}$ is the radalt clearance measurement and $\widehat{x}$ is the state vector that contains positions. A sophisticated search algorithm, designed for the acquisition phase, effectively analyzes topographical data to determine safe flight paths while reducing processing demands.

To assess the efficacy and resilience of the proposed system, we use empirical analysis bolstered by FTI data and comprehensive simulation outcomes. This dual approach highlights the practical application of our AI-enhanced TAN and illustrates its capability to satisfy the stringent certification standards necessary for sophisticated avionic systems. Our study establishes a foundation for a next-generation terrain-aided navigation system that is both safe and reliable by connecting simulation with real-world validation.

While the proposed system builds upon recent advances in adaptive filtering and AI-assisted supervisory control, it is important to acknowledge the foundational contributions in the field. Notably, the Terrain Contour Matching (TERCOM) algorithm [1] and the Sandia Inertial Terrain Aided Navigation (SITAN) system [2,3] have laid the groundwork for batch and recursive TAN architectures, respectively. Our study builds upon these classical methodologies, extending them through hybrid filter banks and artificial intelligence to support robust navigation in GPS-denied environments.

Real-world accuracy of the proposed TAN algorithm to achieve %89 better performance in comparison to just standalone traditional filter solutions. Our extended Kalman filter has better performance than inertial navigation system only solution by %38.

Our contributions in this study are summarized as follows: Development of a hybrid TAN framework and demonstration of reduced navigational drift in environments lacking GPS signals. Our main contribution is the validation of an artificial intelligence-based supervisory system that dynamically selects appropriate navigation algorithms according to terrain characteristics and filter attributes, against real world flight test instrumentation data.

The rest of the paper is organized as follows: Section 2 provides a brief history of previous works. Section 3 explains the TAN methodology. Section 4 introduces our proposed TAN algorithms and Section 5 presents the simulation results. Section 6 concludes our paper.

2. Previous Work

TAN is one of the methods that produces a navigation solution independently of GPS. There are many TAN algorithms in the literature and the most famous ones are terrain contour matching (TERCOM)[2], Sandia inertial terrain aided navigation (SITAN) [3,4] algorithms. TAN algorithms can be examined under three headings: batch, recursive, and hybrid algorithms. Although the first studies on TAN focused on batch or recursive algorithms, current studies generally consist of hybrid algorithms that combine different algorithms, as in this article. Moreover, gravity-aided and magnetic navigation studies are among the alternative methods that can work without GPS.

2.1. Batch TAN Algorithms

In the batch TAN algorithms, the data is processed as a single group. In the TERCOM [2] basically height error comparison is made using different comparison methods. In [5], it is stated that mean absolute difference (MAD) is the most efficient and most correct combination among comparison method for the real time applications. [6] added maximum a posterior (MAP) estimator to TERCOM and SITAN, and noted improvements were achieved.

Probability-based TAN (PTAN) [7] that examines the probability distribution differences is another research. The solution with the minimum probability distribution difference gives the best match. However, [5] stated that the computational load of this algorithm is very high.

[8] observed an increase in system performance using probabilistic data association filter (PDAF) compared to TERCOM. However there are questions marks about running this system in real time [5].

[9] proposes a computer vision based TAN application by using a charge-coupled camera. However, this method has the disadvantage of giving poor results in cases such as adverse weather conditions [5].

2.2. Recursive TAN Algorithms

In the recursive TAN algorithms, the data is processed recursively when new measurement data is received. For SITAN algorithm different filters are compared as nonlinear Kalman filter (KF), extended Kalman filter (EKF), parallel Kalman filter, particle filter (PF). Among these filter structures, PF has the most computational cost. Additionally, it has been observed that a divergence problem may occur in some trials. In this method, terrain linearization is one of the most important factors [4]. In [10] system models have been improved with colored noises for SITAN. Also by adding error models to SITAN, better and more realistic results are achieved [11]. In a study, terrain altitude linearization is made [12]. It is stated that SITAN is used for AFTI/F16 [13] and US Army UH-1 Helicopter [14]. Multi model adaptive estimation (MMAE) and KF is used for this application [14].

In [15] Inertial Terrain Aided Navigation (BITAN)-II method is presented where position and speed errors are examined. Also for navigation solution Bank of KF are used. BITAN [16] which can work in both acquisition and track modes, is actually an improved version of SITAN.

In [17], the point mass filter reached the Cramer-Rao lower limit after convergence, and good results were achieved with this Bayesian method on both flat and variable terrain. However, there are computational problems due to dense meshes [5]. [18] proposed Rao-Backwellisation technique that uses recursive Monte Carlo (MC) runs and particle filter for position estimation.

Viterbi algorithm terrain aided navigation (VATAN) [19] used Viterbi algorithm and contains dynamic programming for position estimate.

Landmark detection is made [20] by using probabilistic data association filter (PDAF). Moreover, a batch implementation is made for PDAF [8].

There are also classical TAN algorithms that perform INS correction using DTED [21,22]. In one study, it was mentioned that it provides a high level of localization accuracy [23]. In addition, there are some image based TAN applications [24,25].

2.3. Hybrid TAN Algorithms

Hybrid TAN algorithms are compose of batch and recursive algorithms. Since it contains recursive algorithms, it is also possible to consider this as recursive algorithms. Hybrid TAN algorithms uses KF and comparison techniques together, thus taking advantage of both batch and recursive methods. With this method it is possible to achieve better navigation solutions [26].

In article [27] written on improved TERCOM, it is mentioned that velocity correction is made as well as position correction using Kalman filter. Via incorporating this method, improvement in accuracy and stability in the navigation solution has been achieved. It has also been stated that this method is a solution to the divergence problem encountered in TERCOM.

TERPROM^® [28], which is used in military systems, produces navigation solutions using INS and radar measurements and digital terrain data. TERPROM^®, which can also operate at low altitudes, has functions such as advanced terrain avoidance cueing (ATAC), predictive ground collision system (PGCAS), and terrain awareness display (TAD).

Recent studies have been proposed to enhance classical TAN systems, including hybrid filters such as batch + APF structures [29], visual TAN using structure-from-motion [30], or deep learning-based object recognition with aerial maps [31]. In particular, marginalized particle filtering (MPF) has become a popular approach to address the nonlinear nature of TAN, with several studies introducing moment-matched kernels [32], boundary flow mechanisms [33], and computationally efficient filter variants [34] to ensure robustness and efficiency. Moreover, adaptively modeling radar altimeter errors [35] and employing LSTM-based measurement models [36] have proven effective in maintaining filter consistency, especially over flat terrain.

High-altitude applications, such as aircraft-based TAN, have also benefited from slant-range measurement strategies using steerable lasers [37], which allow selective targeting of terrain gradients to maximize navigation information. These methods outperform traditional nadir TAN, both in CRLB analysis and simulation, and have been validated with real flight data. Broader perspectives, such as all-source PNT architectures [38,39], highlight the necessity of combining multiple navigation modalities—TAN, SOOP, vision, and sensor fusion—to ensure resilient performance under jamming or spoofing conditions.

TAN is not only used for airborne applications, there are some underwater TAN applications [40,41]. By using sonar and inertial measurement unit (IMU) autonomous navigation of underwater vehicle is investigated. In [42] it is stated that DTED map is not used, terrain height is detected thanks to sonar sensor.

2.4. Gravity-Aided and Magnetic Navigation Algorithms

Apart from TAN, efforts have been made to investigate GPS independent solutions using gravity-aided and magnetism based methods [43,44,45]. Gravity navigation approaches, often implemented in underwater or simulation environments, require precise, often bulky gradiometers and rely on high-resolution gravity maps, making their application on aircraft challenging in terms of both cost and platform dynamics [40,46,47,48]. Moreover, many studies in this domain still focus on improving basic matching or filtering algorithms, indicating that gravity-aided navigation remains in a relatively early stage compared to TAN systems [49]. Similar to this, magnetic navigation techniques rely on anomaly maps and are highly susceptible to magnetic disturbances caused by platforms; this is especially problematic on aircraft because of the complex and dynamic electromagnetic environments present [44,50,51]. While some successful demonstrations have been reported, including real flight tests on fighter aircraft like the F-16 [44], the calibration complexity and variability between platforms make operational deployment difficult.

Biological and marine studies further show how magnetic navigation can be effective in controlled or natural environments, but they also highlight limitations in sensor placement and environmental stability for airborne use [38,52]. However, with off-the-shelf sensors and globally accessible, current terrain databases, TAN systems that use radar altimeters and digital terrain elevation maps provide a more practical option for airplanes. Table 1 provides a comparison of algorithm types.

3. TAN Methodology

3.1. Modeling INS Errors Using Linearized State Equations

In an Inertial Navigation System (INS), the mechanization equations determine the system’s position, velocity, and orientation (attitude) within a specified reference frame. These equations exhibit significant nonlinearity, and sensor flaws, such as bias and noise, eventually degrade navigation accuracy. A typical method to model and mitigate these errors is to linearize the INS error dynamics around a nominal (presumed) trajectory. The resulting set of linear time-varying equations enables the implementation of filtering techniques, including the Kalman filter, for state estimation and error corrections [53].

Given that it is significantly easier for individuals on the Earth’s surface to conceptualize a navigation solution within the local level frame, and considering that the computational errors in the horizontal frame are bounded (Schuler effect), it is deemed more suitable to resolve the error dynamics equations in the local level frame for practical applications.

To formulate the error-state equations, applying a small perturbation to the state variables of the INS in its nominal trajectory, while excluding higher-order terms and negligible terms such as the terms that are multiplied by the square of the Earth’s radius, a system of linear differential equations governing the INS error dynamics can be derived. Furthermore, integrating the sensor error characteristics into this INS error dynamics with specific assumptions, such as first-order Gauss-Markov and random walk, provides the augmented state space equation. Now, using this knowledge, we will derive a linear system model of the errors of INS. Comprehensive knowledge about linearization, INS mechanization and error dynamics derivations are available in [53], Chapters 2, 5 and 6.

The INS mechanization equations can be explained as a set of nonlinear differential equations constructed using IMU measurements taken by a 6-degree-of-freedom system with respect to a specified reference coordinate frame. The general form of the nonlinear ODE can be given as below:

$$\dot{x}\left(t\right)=f(x\left(t\right),u\left(t\right))$$

(1)

Where $x\left(t\right)$ is the INS states and $u\left(t\right)$ is the IMU inputs.

Let the state vector can be composed as

$${\mathbf{x}}_{\mathrm{true}}={\mathbf{x}}_{\mathrm{nominal}}+\delta \mathbf{x},$$

(2)

where ${\mathbf{x}}_{\mathrm{nominal}}$ is the major part and $\delta \mathbf{x}$ is the minor and stochastic part. It can be thought ${\mathbf{x}}_{\mathrm{nominal}}$ as the best guess of the true state (given our motion model and IMU inputs). The model we use for estimation is always subject to some modeling errors and random noises; $\delta \mathbf{x}$ is attributed to the small, stochastic components of the solution.

By taking this into account, linearizing the INS mechanization equations around a nominal solution, we establish that:

$$\delta \dot{x}\left(t\right)={\displaystyle \frac{\partial f\left(x\right(t),u(t\left)\right)}{\partial x}}\delta x+G\left(t\right)={\displaystyle \frac{\partial f\left(x\right(t),u(t\left)\right)}{\partial u}}\delta u$$

(3)

By letting

$$\begin{array}{cc}\hfill F\left(t\right)& ={\displaystyle \frac{\partial f\left[x\right(t\left)\right]}{\partial x}}\mathrm{and}G\left(t\right)={\displaystyle \frac{\partial f\left[x\right(t),u(t\left)\right]}{\partial u}}\hfill \end{array}$$

(4)

we obtain the desired linearized form.

In the local-level frame, the mechanization equations yield a 15-dimensional error-state vector. The position error vector, denoted by $\delta \mathbf{r}={[\delta L,\delta \lambda ,\delta h]}^T$, represents errors in geodetic latitude $\delta L$, longitude $\delta \lambda$ , and altitude $\delta h$. The velocity error vector, $\delta \mathbf{v}={[\delta v_n,\delta v_e,\delta v_d]}^T$, accounts for errors in the north $\delta v_n$, east $\delta v_e$, and down $\delta v_d$ velocity components. The attitude error vector, $\delta \mathit{ϵ}={[\delta r,\delta p,\delta a]}^T$, captures the small-angle errors in roll $\delta r$, pitch $\delta r$, and azimuth $\delta a$. The gyroscope error vector, $\delta \mathit{\omega }={[\delta {\omega }_x,\delta {\omega }_y,\delta {\omega }_z]}^T$, corresponds to gyroscope biases, and the accelerometer error vector, $\delta \mathbf{f}={[\delta f_x,\delta f_y,\delta f_z]}^T$, represents accelerometer biases resolved in x-y-z triad. Concatenating these sub-vectors results in the 15-dimensional error-state vector $\delta {\mathbf{x}}_{15\times 1}^T={\left[\begin{array}{c}\delta \mathbf{r}\delta \mathbf{v}\delta \mathit{ϵ}\delta \mathit{\omega }\delta \mathbf{f}\end{array}\right]}^T.$ After linearizing the INS trajectory along nominal path, the linear relationship between error states and the change of error states can be expressed as follows.

(5)

Where ${\beta }_{\omega x}$, ${\beta }_{\omega y}$,${\beta }_{\omega z}$ and ${\beta }_{fx}$, ${\beta }_{fy}$, ${\beta }_{fz}$ are the autocorrelation time constants of the gyroscope and accelerometer, respectively, resolved in the x-y-z triad., ${\sigma }_{\omega x}$, ${\sigma }_{\omega y}$ and ${\sigma }_{\omega z}$ are the standart deviations associated with gyroscope, ${\sigma }_{fx}$, ${\sigma }_{fy}$ and ${\sigma }_{fz}$ are the standart deviations associated with accelerometer, $R_b^l$ represents the rotation matrix elements which converts body frame to local level frame, $f_n$, $f_e$ and $f_d$ are the accelerometer readings resolved in local NED frame and $w\left(t\right)$ is the unit variance Gaussian white noise.

We can express this set of linear differential equations in an compact way for filtering tasks like state estimation.

$$\dot{\delta \mathbf{x}}\left(t\right)={\mathbf{F}}_x\left(t\right)\delta \mathbf{x}\left(t\right)+\mathbf{G}\mathbf{w}\left(t\right),$$

(6)

${\mathbf{F}}_x\left(t\right)$ is called the dynamic matrix which propagates the INS errors over time, $\mathbf{G}$ is the noise distribution matrix.

3.2. TAN Fundamentals

Terrain-Aided Navigation (TAN) entails determining the platform’s position by comparing the measured terrain profile with the onboard elevation data. These topographical profiles are generated utilizing data from the platform’s barometric altimeter and radar altimeter (RADALT). As depicted in Figure 3, the barometric altimeter provides altitude information relative to mean sea level $\left({\widehat{h}}_{\mathrm{baro}}\right)$, whilst the RADALT measures the platform’s elevation above the ground $\left(h_{\mathrm{radalt}}\right)$. The terrain elevation at the measurement location is determined by subtracting the RADALT value from the barometer altitude.

Batch and recursive solutions have become popular in academia and industry for providing navigation solutions with predictive algorithms using terrain data. The most popular of these algorithms are TERCOM and SITAN, respectively explained in references [2,4]. The SITAN method employs EKF and local terrain linearization algorithms to analyze the terrain profile, with position updates performed as measurements become available. While favorable outcomes may be achieved in certain terrain types and mission profiles, local linearization may not comprehensively address the issue due to the discontinuous and highly nonlinear nature of terrain data in certain areas.

3.3. Terrain Server

Terrain server is a database query system that interpret DTED and provides digital elevation models (DEM). DTED models that is interpreted is a DTED-2 model that has 30m spacing between grid lines. It also has 1 second latitude and longitude spacing. The database for the terrain server is stored onboard computer of the aircraft. It can provide responses with user datagram protocol (UDP) and distributed interactive simulation (DIS). To It has several different functions that handles the LoS, HaT and HoT.

For HaT and HoT queries server uses the provided data to respond with the terrain elevation of a location. It can provide this service for a single location or for a $V\times V$ matrix centered around a location, which is mostly used for our case. LoS is used to predicted the collisions between the aircraft path and the terrain. A straight line is calculated between the aircraft’s location and the stated location, and checked against the elevation of the terrain along the path. To provide these services in a real time environment cache system is used to decrease the time required to load the data into the memory. It stores the DTED in the database to the local memory and retrieves information based on these data.

3.4. Filter Formulations

The formulations of the filters used during implantation are based on the formulations described in this section and are included in the text for completeness.

3.4.1. Kalman Filter Formulation

The Kalman filter is a recursive solution to the linear quadratic estimation problem. It estimates the state ${\mathbf{x}}_k\in {\mathbb{R}}^n$ of a discrete-time controlled process governed by the linear stochastic difference equation

$${\mathbf{x}}_k={\mathbf{F}}_{k-1}{\mathbf{x}}_{k-1}+{\mathbf{G}}_{k-1}{\mathbf{u}}_{k-1}+{\mathbf{w}}_{k-1},$$

(7)

where ${\mathbf{x}}_k$ is the state vector at time step k, ${\mathbf{F}}_{k-1}$ is the state transition matrix, ${\mathbf{G}}_{k-1}$ is the control input matrix, ${\mathbf{u}}_{k-1}$ is the control vector, and ${\mathbf{w}}_{k-1}$ denotes the process noise, assumed to be zero mean Gaussian with covariance ${\mathbf{Q}}_{k-1}$.

The measurement vector ${\mathbf{z}}_k\in {\mathbb{R}}^m$ is defined by

$${\mathbf{z}}_k={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{v}}_k,$$

(8)

where ${\mathbf{H}}_k$ represents the observation matrix and ${\mathbf{v}}_k$ is the measurement noise, also assumed to be zero mean Gaussian with covariance ${\mathbf{R}}_k$.

The Kalman filter comprises two phases: prediction and update.

In the prediction step, the state and error covariance matrix is forecasted using

$${\widehat{\mathbf{x}}}_{k|k-1}={\mathbf{F}}_{k-1}{\widehat{\mathbf{x}}}_{k-1|k-1}+{\mathbf{G}}_{k-1}{\mathbf{u}}_{k-1},$$

(9)

$${\mathbf{P}}_{k|k-1}={\mathbf{F}}_{k-1}{\mathbf{P}}_{k-1|k-1}{\mathbf{F}}_{k-1}^T+{\mathbf{Q}}_{k-1}.$$

(10)

During the update step, the Kalman gain is computed by

$${\mathbf{K}}_k={\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T{\left({\mathbf{H}}_k{\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T+{\mathbf{R}}_k\right)}^{-1},$$

(11)

after which the state estimate is updated with the new measurement ${\mathbf{z}}_k$ according to

$${\widehat{\mathbf{x}}}_{k|k}={\widehat{\mathbf{x}}}_{k|k-1}+{\mathbf{K}}_k\left({\mathbf{z}}_k-{\mathbf{H}}_k{\widehat{\mathbf{x}}}_{k|k-1}\right),$$

(12)

and finally, the error covariance is updated as

$${\mathbf{P}}_{k|k}=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_{k|k-1}.$$

(13)

3.4.2. Unscented Kalman Filter (UKF) Formulation

The Unscented Kalman Filter addresses nonlinear state and measurement models by employing the unscented transform. Let the state vector be ${\mathbf{x}}_k\in {\mathbb{R}}^n$, with a prior estimate ${\widehat{\mathbf{x}}}_{k-1|k-1}$ and covariance ${\mathbf{P}}_{k-1|k-1}$. A set of $2n+1$ sigma points is generated as below, where $\lambda$ is a scaling parameter.

$${\chi }_{k-1}^{\left(0\right)}={\widehat{\mathbf{x}}}_{k-1|k-1},$$

(14)

$${\chi }_{k-1}^{\left(i\right)}={\widehat{\mathbf{x}}}_{k-1|k-1}+{\left(\sqrt{(n+\lambda ){\mathbf{P}}_{k-1|k-1}}\right)}_i,i=1,\dots ,n,$$

(15)

$${\chi }_{k-1}^{\left(i\right)}={\widehat{\mathbf{x}}}_{k-1|k-1}-{\left(\sqrt{(n+\lambda ){\mathbf{P}}_{k-1|k-1}}\right)}_{i-n},i=n+1,\dots ,2n,$$

(16)

In the Prediction step, each sigma point is propagated through the nonlinear state transition function $f(·)$:

$${\chi }_{k|k-1}^{\left(i\right)}=f\left({\chi }_{k-1}^{\left(i\right)},{\mathbf{u}}_{k-1}\right),i=0,\dots ,2n.$$

(17)

The predicted state and covariance are computed as

$${\widehat{\mathbf{x}}}_{k|k-1}=\sum _{i=0}^{2n}{W}_m^{\left(i\right)}{\chi }_{k|k-1}^{\left(i\right)},$$

(18)

$${\mathbf{P}}_{k|k-1}=\sum _{i=0}^{2n}{W}_c^{\left(i\right)}\left({\chi }_{k|k-1}^{\left(i\right)}-{\widehat{\mathbf{x}}}_{k|k-1}\right){\left({\chi }_{k|k-1}^{\left(i\right)}-{\widehat{\mathbf{x}}}_{k|k-1}\right)}^T+{\mathbf{Q}}_{k-1},$$

(19)

where $W_m^{\left(i\right)}$ and $W_c^{\left(i\right)}$ are the weights for the mean and covariance, and ${\mathbf{Q}}_{k-1}$ is the process noise covariance.

In the Update step, each predicted sigma point is transformed through the measurement function $h(·)$:

$${\gamma }_k^{\left(i\right)}=h\left({\chi }_{k|k-1}^{\left(i\right)}\right),i=0,\dots ,2n.$$

(20)

The predicted measurement and its covariance are obtained by

$${\widehat{\mathbf{z}}}_k=\sum _{i=0}^{2n}{W}_m^{\left(i\right)}{\gamma }_k^{\left(i\right)},$$

(21)

$${\mathbf{S}}_k=\sum _{i=0}^{2n}{W}_c^{\left(i\right)}\left({\gamma }_k^{\left(i\right)}-{\widehat{\mathbf{z}}}_k\right){\left({\gamma }_k^{\left(i\right)}-{\widehat{\mathbf{z}}}_k\right)}^T+{\mathbf{R}}_k,$$

(22)

with ${\mathbf{R}}_k$ representing the measurement noise covariance. The cross-covariance between the state and the measurement is computed as

$${\mathbf{P}}_{xz,k}=\sum _{i=0}^{2n}{W}_c^{\left(i\right)}\left({\chi }_{k|k-1}^{\left(i\right)}-{\widehat{\mathbf{x}}}_{k|k-1}\right){\left({\gamma }_k^{\left(i\right)}-{\widehat{\mathbf{z}}}_k\right)}^T.$$

(23)

The Kalman gain is then given by

$${\mathbf{K}}_k={\mathbf{P}}_{xz,k}{\mathbf{S}}_k^{-1}.$$

(24)

Finally, the state estimate and error covariance are updated as

$${\widehat{\mathbf{x}}}_{k|k}={\widehat{\mathbf{x}}}_{k|k-1}+{\mathbf{K}}_k\left({\mathbf{z}}_k-{\widehat{\mathbf{z}}}_k\right),$$

(25)

$${\mathbf{P}}_{k|k}={\mathbf{P}}_{k|k-1}-{\mathbf{K}}_k{\mathbf{S}}_k{\mathbf{K}}_k^T.$$

(26)

3.4.3. Particle Filter (PF) Formulation

The particle filter approximates the posterior density

$$p\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})$$

by a set of N weighted particles

$${\{{\mathbf{x}}_k^{\left(i\right)},w_k^{\left(i\right)}\}}_{i=1}^N,$$

where ${\mathbf{x}}_k^{\left(i\right)}\in {\mathbb{R}}^n$ is a state candidate and $w_k^{\left(i\right)}$ is the corresponding weight.

Assuming a proposal distribution $q\left({\mathbf{x}}_k\right|{\mathbf{x}}_{0:k-1},{\mathbf{z}}_{1:k})$, the weights are updated as

$$w_k^{\left(i\right)}\propto w_{k-1}^{\left(i\right)}{\displaystyle \frac{p\left({\mathbf{z}}_k\right|{\mathbf{x}}_k^{\left(i\right)})p\left({\mathbf{x}}_k^{\left(i\right)}\right|{\mathbf{x}}_{k-1}^{\left(i\right)})}{q\left({\mathbf{x}}_k^{\left(i\right)}\right|{\mathbf{x}}_{0:k-1},{\mathbf{z}}_{1:k})}}.$$

(27)

A common choice is

$$q\left({\mathbf{x}}_k\right|{\mathbf{x}}_{0:k-1},{\mathbf{z}}_{1:k})=p\left({\mathbf{x}}_k\right|{\mathbf{x}}_{k-1}^{\left(i\right)}),$$

(28)

which simplifies the weight update to

$$w_k^{\left(i\right)}\propto w_{k-1}^{\left(i\right)}p\left({\mathbf{z}}_k\right|{\mathbf{x}}_k^{\left(i\right)}).$$

(29)

The weights are normalized by

$$w_k^{\left(i\right)}={\displaystyle \frac{w_k^{\left(i\right)}}{{\sum }_{j=1}^N{w}_k^{\left(j\right)}}}.$$

(30)

The effective sample size is defined as

$$N_{\mathrm{eff}}={\displaystyle \frac{1}{{\sum }_{i=1}^N{\left(w_k^{\left(i\right)}\right)}^2}}.$$

(31)

Resampling is performed when $N_{\mathrm{eff}}<N_{\mathrm{th}}$, where $N_{\mathrm{th}}$ is a predefined threshold, by drawing N particles with replacement according to the weights, and resetting the weights to

$$w_k^{\left(i\right)}={\displaystyle \frac{1}{N}},i=1,\dots ,N.$$

(32)

4. Proposed Algorithms

4.1. Bank of Kalman Filters

In our TAN system, in order to prevent problems caused by high initial errors, different terrain structures affect our altimeter readings and linearization errors, a parallel Kalman filter structure was used , and an appropriate navigation solution was provided by making hard decision making between these filters by artificial intelligence. Here, when the three filter structures that have been explained calculate their separate solutions, the supervisor decides which filter will give better results by looking at certain metrics from the covariance matrices of these solutions, such as trace, determinant, eigenvalue, and the roughness of the terrain profile. The supervisor makes this decision by separating the latitude and longitude solutions of the filters. In other words, in each iteration, supervisor can choose the longitudinal error of one filter while choosing the latitudinal error of the other filter. In this way, by decoupling, errors that may differentially affect the latitude and longitude solutions can be eliminated. Details of EKF, stochastic linearization and parallel filter structures can be examined in [4,5]. Details of UKF and PF are given in [5,54,55].

The terrain data used was converted into a continuous data set using the spatial interpolation method.

4.2. Our TAN Algorithm

We choose the open-loop architecture due to the algorithm’s solution being significantly influenced by the terrain profile, which can produce inaccurate results in steep or non-unique terrains. Moreover, improperly correcting errors in a reliable INS solution, which is already equipped with high-precision, military-grade sensors, creates a risk in flight-critical scenarios. The TAN algorithm (see Figure 4) corrects the errors of INS in an open-loop configuration. It makes positional and velocity update by subtracting the estimated error states from the INS predictions in the iteration TAN runs. In this configuration, unlike a closed loop, the navigation solution at the output has no direct relationship with the next iteration. Bank of filters produce error state estimation outputs in the iteration where the radar measurement is received. The AI supervisor also decides which filter structures output to use. As parallel filter structures, EKF with numeric jacobian, classic unscented Kalman filter and particle filter were preferred.

The primary objective of creating this structure is to minimize the risk of finding incorrect solutions in situations where error states are assumed to follow a Gaussian distribution, thereby allowing the numerical Jacobian retrieval process and the unscented transform to perform well in varying contexts, while simultaneously approximating the distribution of particles affected by noise across various terrain profiles. The measurement of the filter can be expressed as:

$$y_k=h_{radalt}={\widehat{h}}_{baro}-h(\widehat{L_k},\widehat{{\lambda }_k})$$

(33)

In Equation (33), $y_k$ is the estimated terrain height, $h(\widehat{L_k},\widehat{{\lambda }_k)}$ is the estimated height of the terrain obtained from the aircraft’s latitude and longitude using the onboard terrain server database. Our proposed measurement function inserts the adjusted latitude and longitude values into the terrain server to estimate the height of terrain at the true position hypotesis. The input can be formulated as below:

$$\left(\begin{array}{c}\widehat{L}\\ \widehat{\lambda }\end{array}\right)=\left(\begin{array}{c}{L}_{\mathrm{INS}}+\delta {\widehat{L}}_{\mathrm{Filter}}\\ {\lambda }_{\mathrm{INS}}+\delta {\widehat{\lambda }}_{\mathrm{Filter}}\end{array}\right)$$

(34)

In Equation (34), $\left(\begin{array}{c}\widehat{L}\\ \widehat{\lambda }\end{array}\right)$ is the measurement innovation vector between the INS predictions and the estimates of the filter. Where, $\widehat{\delta L}$ and $\widehat{\delta \lambda }$ are the latitude and longitude measurement error estimations of the error state Kalman Filter.

4.3. AI Supervisor

One way to have more precise navigation is to have better sensors or better algorithms; however, if sensors and algorithms are fixed beforehand, one can cleverly choose which algorithm serves better in which terrain conditions. This is how our AI supervisor is conceived.

As it can be observed from Figure 4, our AI supervisor takes two series of inputs: the covariance matrices of bank of filters and the roughness values of the terrain. Supervisor is a classifier model which selects the predicted best rms value returning filter.

The output of the model are two label structures consisting of three classes separately on the latitude and longitude axes. Thus, a total of nine different filter configurations are obtained, and the model selects the most suitable one from these nine.

To select our model architecture, we used hyper parameter optimization tool Optuna. Layer count, dropout layer batch size, hidden state size are all given as a parameter to our optimization pipeline. Our final model has 3 linear layers to process the input. Each layer has a hidden with fully connected (FC) layers size of 1024 and has RELU non linear activation function. The architecture of the supervisor model is given in Figure 5.

13 different terrains were selected, including rough, moderate and flat terrain. 10 Monte Carlo flights of 1000 seconds were performed on these terrains. Filter error metrics, covariance data and roughness data where the flight took place were recorded. The model was then trained on this dataset collected from simulations in the environment mentioned in the Section 5.1. 85% of the data was used as training data and 15% was used as test data. Training was carried out for 30 epochs using a cross-entropy loss function applied to both outputs. The model achieved a validation accuracy of 0.88, computed over the joint label space consisting of all nine class combinations. The Figure 6 illustrates the learning dynamics of the proposed architecture.

5. Simulations and Results

5.1. Simulations

In this section, simulation environment designed to provide AI training by generating the results of the aforementioned methods and filters will be explained. First of all, three different terrain models were selected for simulations: rough, moderate and flat terrain. Then, flight trajectories were determined. While determining these trajectories path height standard deviation ${\sigma }_T$ and path height slope standard deviation ${\sigma }_Z$ parameters are examined. Even though both parameters are a metric of roughness, ${\sigma }_Z$ is more significant metric in such filtering methods since algorithms such as KF use linearization methods [5]. ${\sigma }_T$ corresponds to the standard deviation of the terrain heights on the trajectory. N is the samples, $h_i$ is the terrain height measurement and $\overline{h}$ is the mean height of that terrain.

$$\overline{h}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{h}_i$$

(35)

By using $h_i$ and $\overline{h}$ path height standard deviation ${\sigma }_T$ can be calculated as below.

$${\sigma }_T=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{(h_i-\overline{h})}^2}$$

(36)

In order to calculate the path height slope standard deviation ${\sigma }_Z$, slopes $d_i$ should be found first and $\overline{d}$ is the mean of the slopes.

$$d_{i-1}=h_{i-1}-h_i$$

(37)

$$\overline{d}={\displaystyle \frac{1}{N-1}}\sum _{i=1}^{N-1}{d}_i$$

(38)

By using $d_i$ and $\overline{d}$ path height slope standard deviation ${\sigma }_Z$ can be calculated as below.

$${\sigma }_Z=\sqrt{{\displaystyle \frac{1}{N-1}}{\sum }_{i=1}^{N-1}{(d_i-\overline{d})}^2}$$

(39)

Sample ${\sigma }_T$ and ${\sigma }_Z$ values of chosen terrains are given in the Table 2.

After the flight trajectories were selected, 1000-second flights were made on the determined paths. A trajectory sample is given in Figure 7. A path was defined for the aircraft such that there was a random acceleration on the horizontal axis. Then, the IMU was modeled to make measurements. The characteristics of the accelerometer, gyroscope and barometric altimeter used in the simulations are given in Table 3.

5.2. Validation with Real World FTI Data

The real flight test data (FTI) of Hürjet advanced jet trainer platform projects within Turkish Aerospace was used. The TAN algorithm was run at a frequency of 20 Hz. The data were recorded with a sampling frequency of 128 Hz. There are three different navigation solutions: free inertial solution (FI), GPS solution, and hybrid solution. The hybrid solution is a filtered version of the FI and GPS solutions and produces more accurate and reliable results compared to both solutions. Therefore, hybrid solution was accepted as "ground truth" while evaluating the performance of TAN algorithm. FI navigation solution shows drift over time due to various errors as explained in Section 3.1 and constitutes the input data that are being tried to improve. Similarly, radalt measurements were also taken from the aircraft. When the aircraft makes high angles of roll or yaw, radalt cannot measure the HaT correctly and produces invalid data. For this reason, when the radalt measurement could not be produced, only the prediction step was applied and no correction was made.

5.3. Results

In order to compare to results of the trials two parameters are examined: absolute error and circular error probable (CEP). Absolute error and CEP [56] are computed as below.

$${\mathbf{Error}}_i=\sqrt{e_{x_i}^2+e_{y_i}^2}$$

(40)

$${\mathbf{C}}_0=\left\{\begin{array}{cc}0.562{\sigma }_x+0.615{\sigma }_y,\hfill & {\sigma }_y<{\sigma }_x\mathrm{and}{\displaystyle \frac{{\sigma }_y}{{\sigma }_x}}<0.28\hfill \\ 0.562{\sigma }_y+0.615{\sigma }_x,\hfill & {\sigma }_x<{\sigma }_y\mathrm{and}{\displaystyle \frac{{\sigma }_x}{{\sigma }_y}}<0.28\hfill \\ 1.1774{\displaystyle \frac{({\sigma }_x+{\sigma }_y)}{2}},\hfill & \mathrm{else}\hfill \end{array}\right.$$

(41)

$$CEP=\sqrt{C_0^2+{\mu }_x^2+{\mu }_y^2}$$

(42)

Here $e_{x_i}$ and $e_{y_i}$ are the position errors between the ground truth and the TAN output. The standard deviations of the position errors between the ground truth and the TAN output is denoted as ${\sigma }_x$ and ${\sigma }_y$. Also ${\mu }_x$ and ${\mu }_y$ are the means of the position errors between the ground truth and the TAN output. CEP is calculated by using coefficient $C_0$.

The FTI data absolute error results are given in Figure 8 and Figure 9 for two different flights. Parallel filter solutions give better navigation solutions than INS only solutions as proposed.

In Figure 8 and Figure 9 The reason for the dips in the Estimated Line between 250-300, 480-540, 720-760 and 930-1000 seconds is that the aircraft is performing an agile maneuver at these moments. This can be called radalt shadowing due to roll and pitch manuevers. Moreover, there is a variable called "radalt attitude angle" that is dependent on the roll and pitch angles of the aircraft. When this variable is below a certain value, measurements from radalt can be accepted as correct. On the specified times that listed above, the attitude angle is above threshold and the correction to the bank of filters can not be accepted. In order to visualize this type of manuevers, estimated HaT value is set to zeros. This situations can also be validated from the dip to 0 in the Radalt measurements in the intervals mentioned above.

6. Conclusions

Designed for the Hürjet advanced jet trainer, the paper presents a new artificial intelligence-enhanced Terrain-Aided Navigation (TAN) system. The TAN system improves navigation precision and computation efficiency by including adaptive filtering techniques, an advanced terrain server, and an artificial intelligence supervisory module that dynamically alternately between multiple filtering schemes (EKF, UKF, PF). Comparatively to conventional navigation methods (INS-only solutions), the proposed hybrid TAN method significantly increases low-altitude navigation efficacy. Emphasizing improved navigational accuracy even in GPS-denied environments, Monte Carlo simulations and empirical flight test instrumentation (FTI) data support the efficacy and viability of the AI-driven methodology. With an accuracy of around 82% in appropriate filter selection, the AI supervisor efficiently finds acceptable filters based on terrain roughness and ambient factors, hence significantly lowering absolute error values. This study provides important new insights on the issues and operational aspects related to avionics certification by using real flight data from the Hürjet platform to link theoretical modeling with pragmatic application. We understand the importance of upcoming projects and stress the change from simulation to real-time onboard applications as well as a complete certification process for the TAN integration with aircraft control systems. Appropriate for inclusion into advanced avionics systems, the work represents a notable advancement in autonomous, trustworthy, and accurate terrain-aided navigation procedures.