Analysis of the Trusted Inertial Terrain-Aided Navigation Measurement Function

3.1 Derivation of the MRPD Model

This derivation is based on the work of Johnstone and Shene (1992) on the intersection of a cone and plane. Two approximations are invoked so that the MRPD optimization problem can be formulated as a cone–plane intersection problem with analytic solutions: (1) each MRPD Doppler bin is approximated by a single Doppler value (Doppler cone approximation), and (2) the terrain near each MRPD ground point is approximated by a plane (planar terrain approximation). Note that there is a separate planar terrain approximation for each individual MRPD ground point. The planar terrain approximation is generally valid over distances on the order of a dozen meters, and the Doppler cone approximation is reasonable provided that the width of a Doppler bin (typically several meters) is small relative to the radar’s altitude and to the extent of validity of the planar terrain approximation. It will be shown that the closed-form solution to the MRPD optimization problem lies at the vertex of a hyperbola formed by the intersection of the Doppler cone and planar terrain and that this closed-form solution may be expressed in terms of the navigation states, radar parameters, and terrain structure around the MRPD ground points.

Unless stated otherwise, all vectors are expressed in a local frame ℓ; thus, the superscript is dropped hereafter for convenience of notation.

Let us assume that the MRPD optimization problem has been solved once, that the MRPD ground points have been located, and that local planar approximations have been formed around each MRPD ground point. Let the Doppler cone for a single Doppler bin be parameterized by its vertex p_r ∈ ℝ³ (radar position), a unit vector ${\hat{v}}_r\in S^2$ along the cone’s axis (unit vector of radar velocity), and an angle α ∈ (0, π) from the cone’s axis to the edge of the cone (the Doppler angle approximating the Doppler bin under consideration). Furthermore, let us suppose that the terrain plane near the associated MRPD ground point is parameterized by p_g ∈ ℝ³ (the solved MRPD ground point) and a unit vector $\hat{o}\in S^2$ normal to the plane (the terrain normal at the solved MRPD ground point). Let t ∈ ℝ³ be the intersection point between the cone’s axis and the plane (assuming it exists), and let d ∈ ℝ³ be the vector from p_r to t parameterized by the scalar t ∈ ℝ as follows:

$d=t-p_r=t{\hat{v}}_r$4

Following Johnstone and Shene (1992), we use the following definitions:

${\hat{v}}_r^{\text{'}}=\{\begin{array}{cc}{\hat{v}}_r& \text{if}t\ge 0\\ -{\hat{v}}_r& \text{otherwise}\end{array}$ 5

and:

${\hat{o}}^{\prime }=\{\begin{array}{cc}\hat{o}& \text{if}{\hat{v}}_r^{\prime T}\hat{o}\ge 0\\ -\hat{o} & \text{otherwise}\end{array}$6

so that ${\hat{v}}_r^{\prime }$ points from the radar to the intersection point t and the angle θ between ${\hat{v}}_r^{\prime }$ and ${\hat{o}}^{\prime }$ is always positive and acute: θ ∈ (0, π/2). The parameters θ, ${\hat{v}}_r^{\prime }$, and ${\hat{o}}^{\prime }$ are then related as follows:

$\cos \left[\theta \right]={\hat{v}}_r^{\prime T}{\hat{o}}^{\prime }$7

The intersection of a plane and an out-of-plane cone forms either a parabola, hyperbola, ellipse, or circle, depending on the angles θ and α. Of these four cases, the hyperbola case is of general interest. The circle and ellipse cases manifest when cos[θ] > sin[α], the parabola case manifests when cos[θ] = sin[α], and the hyperbola case manifests when cos[θ] < sin[α]. The circle and ellipse cases generally occur with small Doppler angles (α ≈ 0 or α ≈ π) or when the radar is plummeting toward the ground (θ ≈ 0). The parabola case only occurs under a singular geometric condition with probability zero. If we limit our consideration to cases in which the radar is flying with a moderate pitch angle relative to the ground, the radar is pointed towards the nadir, and the radar beam width is moderate (< 45°), the condition cos[θ] < sin[α] is often satisfied and only the hyperbola case needs be solved. In these cases, the MRPD optimization problem is reduced to minimizing the distance between the cone’s center and a point along the resulting hyperbola. Figure 2 illustrates the geometry for the hyperbola case.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Cone-plane intersection geometry for p_r = [0, 0, 15]^T, ${\hat{v}}_r={[1,0,0]}^T,\hat{o}={[0,0,1]}^T$, p_g = [1, 1, 0]^T, and α = 45°. The MRPD point lies at the vertex of the Doppler bin hyperbola in front of the radar. An aliased hyperbola resulting from the two-way definition of the cone lies behind the radar, although it is of no interest here.

Although it is certainly possible for one of the other cases to manifest, such cases are rare when the flight pitch angle is moderate and the radar is pointed towards the nadir. In fact, in the Monte Carlo simulations in the following sections, no ellipse or parabola cases manifested. Still, depending on the situation, the ellipse case may need to be solved; this case is deferred to future work.

There are two minimizers for this hyperbolic distance minimization problem, one at each hyperbola vertex (a corresponding proof is given in Appendix A). Because the hyperbola lies in the ground plane, it is natural to parameterize the hyperbola and its associated components with vectors in the plane. Let the projection of ${\hat{v}}_r^{\prime }$ onto the terrain plane be expressed as follows:

${\hat{m}}^{\prime }=\frac{{\hat{o}}^{\prime }\times \left({\hat{v}}_r^{\prime }\times {\hat{o}}^{\prime }\right)}{\Vert {\hat{o}}^{\prime }\times \left({\hat{v}}_r^{\prime }\times {\hat{o}}^{\prime }\right)\Vert }=\frac{{\hat{v}}_r^{\prime }-\left({\hat{v}}_r^{\prime T}{\hat{o}}^{\prime }\right){\hat{o}}^{\prime }}{\Vert {\hat{v}}_r^{\prime }-\left({\hat{v}}_r^{\prime T}{\hat{o}}^{\prime }\right){\hat{o}}^{\prime }\Vert }$8

and let the following be an orthogonal unit vector in the plane such that $\{{\hat{m}}^{\prime },{\hat{n}}^{\prime }$, and ${\hat{o}}^{\prime }\}$ form a right-handed, orthonormal triad:

${\hat{n}}^{\prime }=\frac{{\hat{o}}^{\prime }\times {\hat{v}}_r^{\prime }}{\Vert {\hat{o}}^{\prime }\times {\hat{v}}_r^{\prime }\Vert }$9

Then, a point p_ϕ ∈ ℝ³ on the hyperbola is given by the following:

$p_{\varphi }=c_0\pm \cosh \left[\varphi \right]m{\hat{m}}^{\prime }+\sinh \left[\varphi \right]n{\hat{n}}^{\prime }$10

where c₀ ∈ ℝ³ is the center of the hyperbola (the point on the plane halfway between the two hyperbola vertices), ϕ ∈ ℝ is the hyperbolic parameter, and ± indicates that there are two sides to the hyperbola — one side is of interest, and the other is alias owing to the two-sided definition of the cone (see Figure 2). The two sides may be distinguished by comparing the angle between p_ϕ – p_r and ${\hat{v}}_r$ against α. The quantities $m{\hat{m}}^{\prime }$ and $n{\hat{n}}^{\prime }$ are the semi-major and semi-minor axes of the hyperbola, respectively, where ${\hat{m}}^{\prime }$ and ${\hat{n}}^{\prime }$ are the directions and m ∈ ℝ and n ∈ ℝ are the magnitudes.

The two hyperbola vertices are located at ϕ = 0 and are given as follows:

$p_0=c_0\pm m{\hat{m}}^{\prime }$11

From the work of Johnstone and Shene (1992), the center of the hyperbola is given as follows:

$c_0=p_r+\Vert d\Vert \left({\hat{v}}_r^{\prime }+\frac{\sin \left[\alpha \right]}{2}\left(\frac{1}{\cos [\alpha +\theta ]}-\frac{1}{\cos [\alpha -\theta ]}\right){\hat{m}}^{\prime }\right)$12

with the following semi-major axis length:

$m=-\frac{\Vert d\Vert \sin \left[\alpha \right]}{2}\left(\frac{1}{\cos [\alpha +\theta ]}+\frac{1}{\cos [\alpha -\theta ]}\right)$13

Substituting Equations (12) and (13) into Equation (11) produces the following expression for the hyperbola vertices:

$p_0=p_r+\Vert d\Vert \left({\hat{v}}_r^{\prime }\mp \frac{\sin \left[\alpha \right]}{\cos [\alpha \mp \theta ]}{\hat{m}}^{\prime }\right)$14

The range from the cone’s center to the two vertices is then written as follows:

$r_0=\Vert p_0-p_r\Vert =\Vert d\Vert \Vert {\hat{v}}_r^{\prime }\mp \frac{\sin \left[\alpha \right]}{\cos [\alpha \mp \theta ]}{\hat{m}}^{\prime }\Vert$15

This range equation is simplified in two ways. First, the line-plane distance ||d|| is replaced by the following formula:

$\Vert d\Vert =\left|\frac{{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }}{{\hat{v}}_r^{\prime T}{\hat{o}}^{\prime }}\right|=\left|\frac{{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }}{\cos \left[\theta \right]}\right|$16

Second, ${\hat{v}}_r^{\prime }$ is decomposed into $\left\{{\hat{m}}^{\prime },{\hat{n}}^{\prime },{\hat{o}}^{\prime }\right\}$ components, noting that ${\hat{v}}_r^{\prime }$ and ${\hat{n}}^{\prime }$ are orthogonal:

${\hat{v}}_r^{\prime }=\left({\hat{v}}_r^{\prime T}{\hat{m}}^{\prime }\right){\hat{m}}^{\prime }+\left({\hat{v}}_r^{\prime T}{\hat{o}}^{\prime }\right){\hat{o}}^{\prime }=\sin \left[\theta \right]{\hat{m}}^{\prime }+\cos \left[\theta \right]{\hat{o}}^{\prime }$17

Substituting Equations (16) and (17) into Equation (15) and simplifying produces the following expression for the range to the two hyperbola vertices:

$r_0=\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right| \left|\sec \right[\alpha \mp \theta \left]\right|$18

To summarize, this analytic expression for the minimum range in a prescribed Doppler bin is a function of the radar position p_r, the solved MRPD ground point p_g, the terrain normal at the MRPD ground point $\hat{o}$, the prescribed Doppler angle α, and the angle between the terrain normal and velocity vector θ. Of course, the accuracy of this expression depends on whether the solved MRPD ground point p_g is a reasonable proxy for the actual point of reflection and whether the the local planar terrain model around p_g is an acceptable approximation. Errors arising from these assumptions on p_g and ${\hat{o}}^{\prime }$ are explored in Section 5.

One assumption invoked in this derivation should be addressed: the assumption that the Doppler cone axis intersects the terrain plane. For the case in which no intersection exists (the radar is flying perfectly parallel to the terrain), Johnstone and Shene (1992) provided a separate solution for the resulting hyperbola (Theorem 4.2). Fortunately, the resulting range-to-hyperbola-vertex expression can be recovered from Equation (18) in the limit as θ → π/2.

3.2 Limits of the MRPD Equation

Two limits of the analytical MRPD range model in Equation (18) are explored to verify that the model is reasonable. First, in the limit as θ → π/2, one can recover the Pythagorean distance expression for an aircraft flying parallel to the ground plane at a ground clearance h > 0 and measuring the distance to a point at an angle α ∈ (0, π) from the horizontal — a reassuring check:

$\underset{\theta \to \pi /2}{\lim }{r}_0=\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right|\left|\csc \right[\alpha \left]\right|=\frac{h}{\sin \left[\alpha \right]}$19

Second, in the limit as both θ → π/2 and α → π/2 (the radar is flying parallel to the ground plane and looking straight down), the minimum range reduces to just the ground clearance h — another reassuring check:

$\underset{\theta ,\alpha \to \pi /2}{\lim }{r}_0=\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right|=h$20

3.3 Jacobians of the MRPD Equation

It is straightforward to compute the position Jacobian of Equation (18), as only the first term varies with it:

$\frac{\partial r_0}{\partial p_r}=-\mathrm{sgn}\left[{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right]\left|\sec \right[\alpha \mp \theta \left]\right|{\hat{o}}^{\prime }$21

Here, sgn is the sign function. The velocity Jacobian is more complex, as both α and θ depend on v_r. The velocity Jacobian is as follows:

$\frac{\partial r_0}{\partial v_r}=\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right|\left|\sec \right[\alpha \mp \theta \left]\right|\tan [\alpha \mp \theta ]\left(\frac{\partial \alpha }{\partial v_r}\mp \frac{\partial \theta }{\partial v_r}\right)$22

The partial derivative of the Doppler angle α w.r.t. the velocity vector is determined via implicit differentiation of its defining equation:

$\dot{r}=-\Vert v_r\Vert \cos \left[\alpha \right]$23

Here, $\dot{r}$ ∈ ℝ is the prescribed range rate — a constant — which is related to the radar Doppler shift in Equation (2) by $\Delta f=-2\frac{f_c}{c}\dot{r}$. Differentiating Equation (23) and solving for $\partial \alpha /\partial v_r$ produces the following:

$\frac{\partial \alpha }{\partial v_r}=\frac{\cot \left[\alpha \right]}{\Vert v_r\Vert }{\hat{v}}_r=\mathrm{sgn}\left[t\right]\frac{\cot \left[\alpha \right]}{\Vert v_r\Vert }{\hat{v}}_r^{\prime }$24

Similarly, the partial derivative of the angle θ is determined via implicit differentiation of its defining expression in Equation (7):

$\frac{\partial \theta }{\partial v_r}=-\mathrm{sgn}\left[t\right]\frac{\csc \left[\theta \right]}{\Vert v_r\Vert }\left(I_{3\times 3}-{\hat{v}}_r^{\prime }{\hat{v}}_r^{\prime T}\right){\hat{o}}^{\prime }=-\mathrm{sgn}\left[t\right]\frac{\csc \left[\theta \right]}{\Vert v_r\Vert }\left({\hat{o}}^{\prime }-\cos \left[\theta \right]{\hat{v}}_r^{\prime }\right)$25

Substituting Equations (24) and (25) into Equation (22) and simplifying produces the final expression for the velocity Jacobian:

$\frac{\partial r_0}{\partial v_r}=\mathrm{sgn}\left[t\right]\frac{\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right|\left|\sec \right[\alpha \mp \theta \left]\right|\tan [\alpha \mp \theta ]}{\underset{\left(\left(\cot \right[\alpha ]\mp \cot [\theta \left]\right){\hat{v}}_r^{\prime }\pm \csc \left[\theta \right]{\hat{o}}^{\prime }\right)}{\Vert v_r\Vert }}$26

Equations (21) and (26) represent the first-order variation of the MRPD measurement function with small perturbations in the radar position and velocity, respectively. These equations elucidate two remarkable properties of the measurement function. First, the direction of the position Jacobian is exclusively determined by the terrain normal at the MRPD ground point ${\hat{o}}^{\prime }$, whereas the magnitude is exclusively determined by the angles α and θ. Second, the direction of the velocity Jacobian is determined by a superposition of the terrain normal at the MRPD ground point ${\hat{o}}^{\prime }$ and the velocity unit vector ${\hat{v}}^{\prime }$, and the magnitude of the Jacobian is determined by the angles α and θ, as well as the radar’s speed v_r and approximate ground clearance $\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right|$. Specifically, the magnitude of the velocity Jacobian varies proportionally with the altitude and inversely with the speed: the MRPD measurements become more sensitive to radar velocity perturbations as the radar’s ground clearance increases, and the measurements become less sensitive to the same perturbations as the radar’s speed increases.

This relationship among altitude, speed, and range measurement sensitivity to velocity errors is succinctly captured in Equation (24), the partial derivative of the Doppler angle w.r.t. radar velocity. This equation represents the first-order variation of the Doppler angle with small perturbations in the radar’s velocity. The Doppler angle is important, as its projection onto the ground defines a Doppler bin and, thus, the possible locations of an MRPD ground point. Specifically, at an altitude h > 0, a velocity perturbation Δv ∈ ℝ³ creates an along-track shift Δx ∈ ℝ³of a Doppler bin by approximately $\Delta x\approx h\tan \left[\frac{\cot \left[\alpha \right]}{\Vert v_r\Vert }{\hat{v}}_r^T\Delta v\right]$. This sensitivity of the Doppler angle and Doppler bin location to velocity perturbations will be important in Section 5.

For a “nominal” set of flight conditions θ ∈ [60°, 90°) (less than 30° flight pitch angle relative to the ground) and α ∈ [60°, 120°] (less than 30° Doppler angle from the zero-Doppler angle), the magnitude of the position Jacobian is bounded to the interval [1, 2]. For the same limits, the magnitude of the velocity Jacobian is approximately bounded to the interval [0, 4h/||v_r||], where $h=\left|{\left(p_g-p_r\right)}^T{\hat{o}}^{\prime }\right|$ is the approximate ground clearance. In these cases, the velocity Jacobian direction is primarily aligned with the ${\hat{o}}^{\prime }$ direction as cot[x] ≈ 0 and csc[x] ≈ 1 for x ≈ 90°.

5.1 Synthetic Terrain Generation

The FIM model in Equations (31) and (33)–(35) provides insight into the information extractable from a single MRPD measurement set; however, as explained above, this model does not characterize the general performance of the system. Instead, according to Equation (38), one must average over realizations of Equation (31) to characterize the modeled general performance. This section invokes Monte Carlo simulations to generate realizations of synthetic terrain from which sequences of MRPD measurement sets are extracted. The single-epoch FIMs from Equation (31) associated with these individual MRPD measurement sets are averaged to characterize the general performance of the modeled system. Synthetic terrain is invoked so that individual radar–terrain geometries may be randomized while controlling the overall statistics of the geometry.

Five synthetic terrain maps with different levels of terrain roughness are considered. It is recognized that synthetic terrain cannot replace real terrain for the purposes of quantifying the performance of TAN algorithms. Yet, it is infeasible to evaluate every potential measurement scenario. The five synthetic terrain scenarios enable a controlled verification and validation of the analytical range, Jacobian, and FIM models developed in the previous sections, as well as a characterization of the general behavior of the TITAN algorithm as it responds to terrains of varying roughness under the acknowledged caveat that the evaluated behavior will never perfectly match reality.

There exists little literature on synthesizing terrain for the purposes of evaluating TAN algorithms — certainly no standard exists. Most authors, such as Bergman (1997), do not synthesize terrain and instead evaluate the performance of TAN algorithms over digitized versions of real terrain. Some older TERCOM studies synthesized two-dimensional (2D) Gaussian terrain with a circular autocorrelation function of the form $\exp \left[-\sqrt{x^2+y^2}/{\sigma }_c\right]$, where x ∈ ℝ and y ∈ ℝ are the two horizontal lags and σ_c > 0 is the correlation length (Cannon, 1974; Cannon Jr, 1978). These authors note that the performance of TERCOM over real and synthetic terrain is “almost equivalent” when the statistics of synthetic and real terrain are similar, partially addressing concerns related to the verification of TAN algorithms with synthetic terrain.

A method similar to that of Cannon (1974) was invoked to generate isotropic random surfaces for the terrain of each of the five maps. Specifically, the terrain was generated by convolving 2D zero-mean Gaussian random noise with a circular 2D Gaussian kernel — a method recommended by Garcia and Stoll (1984) for controlled modeling of rough random surfaces. Following the method of Garcia and Stoll (1984), a square DEM of size M × M with side length L > 0, post spacing Δl > 0, elevation standard deviation σ_h > 0, and terrain correlation length σ_c > 0 was generated via the following three steps:

1. Fill a square patch G ∈ ℝ^M×M with independent samples taken from a scalar zero-mean Gaussian distribution with standard deviation σ_h.

2. Create a 2D kernel K ∈ ℝ^M×M of the form $K[x,y]=\frac{2L}{M\sqrt{\pi }{\sigma }_c}\exp \left[-\frac{1}{2}\frac{x^2+y^2}{{\left({\sigma }_c/2\right)}^2}\right]$, where $x,y\in \left[-\frac{L}{2},+\frac{L}{2}\right]$ and $\frac{2L}{M\sqrt{\pi }{\sigma }_c}$ is a scale factor that shapes the statistics of the resulting map.

3. Convolve the random patch and kernel with a normalized Fourier transform to create the isotropic random surface: $S={ℱ}^{-1}\left[ℱ\right[G]\odot ℱ[K\left]\right]$, where ⊙ indicates the Hadamard product and ℱ and ${ℱ}^{-1}$ are the forward and inverse Fourier transforms.

MATLAB code for the synthetic terrain generation is based on that provided by Bergström (2012). The resulting random surface S ∈ R^M×M has several convenient properties. First, the terrain elevation values are normally distributed with zero mean and a variance of ${\sigma }_h^2$. Second, the marginal terrain gradients in the row and column directions are normally distributed with zero mean and a variance of $2{\sigma }_h^2/{\sigma }_c^2$. Third, the terrain normal vectors are symmetrically distributed about the vertical. Finally, these statistics do not spatially vary on wide scales over the map.

Five square, 20-km-wide, synthetic DEMs were generated with this method. Each map was generated using 10,000 posts (uniformly spaced sample points), a post spacing of 2 m, a terrain elevation standard deviation of 10 m, and a terrain correlation length of {50, 100, 150, 200, 250} m. The mean elevation of each map was adjusted such all elevation values were positive. Figure 3 depicts two patches of terrain synthesized with 50- and 250-m correlation lengths. Table 1 provides some statistics related to the terrain elevation and slopes over the five maps. Figure 4 plots the distribution of the angles of the terrain normals away from the vertical to evidence their azimuthally symmetric distribution.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Synthetic isotropic terrain with an elevation standard deviation of σ_h = 10 m and a correlation length of σ_c = 50 m (left) or σ_c = 250 m (right)

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Angular distribution of the terrain normals for the 50-m correlation length map The pitch angle of the terrain normals is uncorrelated with the azimuth angle; thus, the terrain normals have no preference in the azimuthal direction.

The purpose of these terrain maps is to evaluate the range, Jacobian, and CRLB models over a variety of radar–terrain geometries while controlling the statistics of the elements that govern the two expressions. Accordingly, the five terrain correlation lengths were chosen to exercise a variety of θ angles, upon which the analytical expressions depend. As a reference for these five synthetic maps, the median terrain slope on the moon — one proposed application domain of the TITAN algorithm (Haydon & Humphreys, 2023) — is 2°–7°, depending on the specific location on the moon (Rosenburg et al., 2011). Referring to Table 1, this median slope approximately corresponds with the 150-, 200-, and 250-m correlation maps, indicating that these maps might be considered “average” terrain, whereas the 50- and 100-m correlation maps might be considered “extreme” terrain.

5.2 Monte Carlo Setup

The following subsections verify the derived MRPD, Jacobian, and CRLB models using data generated by Monte Carlo simulations. A radar was simulated as flying horizontally at an altitude of 9000 ft (2743 m) with a speed of 80 m/s — a replication of the experimental flight conditions reported by Haydon and Humphreys (2023). The radar was configured with a circular beam width diameter of 25° and a ground-projected Doppler bin width of ≈4 m. The ground-projected radar beam spans ≈1200 m. The Monte Carlo simulation swept over two parameter sets: (1) the terrain correlation lengths σ_c ∈ {50, 100, 150, 200, 250} m and (2) the number of prescribed Doppler bins N ∈ {6, 8, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100}. The Doppler bins were configured so that the corresponding prescribed Doppler angles were uniformly distributed over [–10°, +10°] from the nadir (e.g., three Doppler bins would correspond with the three Doppler angles α ∈ {100°, 90°, 80°}).

The radar’s horizontal position was randomized 500 times for each combination of map and number of Doppler bins, ensuring that the ground-projected beam widths were always contained within the limits of the map. For each instance, a search was performed to find the MRPD ground points and associated range values. During the searches, the elevation maps were linearly interpolated to a 1-m resolution in order to reduce the impact of the discrete nature of the map on the solved range values. For each MRPD ground point, the local terrain normal was estimated with a least-squares solver using all of the terrain patch points in a 10-m-square neighborhood. Finally, the MRPD search was repeated 12 more times while perturbing the radar position by 2 m and velocity by 10 cm/s, and a finite-difference Jacobian of the measurement function was approximated. The MRPD values, ground points, terrain normals, and Jacobians were all recorded for each combination of map and number of Doppler bins for the following analysis.

5.3 Range Model Verification

Equation (18) provides an analytical expression for the MRPD value under thin Doppler cone and planar terrain assumptions. To verify the validity of this expression and, by proxy, these two assumptions, the analytical range expression was compared against the empirically found minimum range values from the Monte Carlo simulations.

Figure 5 presents the cumulative distribution function (CDF) of the difference between analytical MRPD values and the corresponding empirically found MRPD values for each of the five terrain maps. Over all maps, the minimum difference was –340 m, the maximum difference was 2.92 m, and the median difference was approximately 0.14 m. Overall, the empirically found and analytical MRPD values are very close, although there are some outliers. Excepting the 50-m correlation length case, 99.9% of the range differences are smaller than ≈ 1 m. The remaining 0.1% could be attributed to inexact planar approximations or the 1-m discretization of the terrain maps (errors in p_g or $\hat{o}$).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Cumulative distribution of the difference between analytical MRPD values and empirical MRPD values

The CDF is split across the five terrain maps to illustrate how the simulated terrain roughness impacts the accuracy of the analytical range equation. The magnitude of the range differences tends to increase with increasing terrain roughness (decreasing terrain correlation length). As the terrain becomes rougher, the elevation gradients are magnified, the terrain normal vectors deviate more significantly from the vertical, and the local curvature of the terrain increases. These three trends degrade the validity of the planar terrain approximation and the subsequent accuracy of Equation (18). It is unsurprising, then, that the range produced by Equation (18) diverges somewhat from the empirically determined range as the terrain roughness increases.

Several other properties of Figure 5 should be noted. First, there is a long negative tail for the range differences, but only a short positive tail. That is, the analytical expression for the MRPD value sometimes produces a smaller range than that produced by the empirical search, but only rarely produces a larger MRPD value. This asymmetry in the tails can be explained by the fact that the planar terrain approximation is formed using an already determined MRPD ground point. Because the analytical MRPD expression uses this point (p_g), by definition, it already has an upper bound on the MRPD value equal to the range to this point. The opposite is not true — the analytical expression could always yield a smaller range value. The fact that the analytical expression can produce a longer range despite being furnished an alleged upper bound can be attributed to the difference in Doppler bin models between the analytical and empirical approaches. The analytical Doppler bin constraint models the Doppler cone as an infinitely thin cone centered on a real Doppler bin of some width. As a result of this modeling difference, it is possible for the empirical search to select a point in the Doppler bin that is closer to the radar than any point on the analytical cone, yet not technically along the analytical Doppler cone. Stated otherwise, the domains of the two Doppler bins are not the same; hence, it is possible for the two to arrive at different solutions.

Second, the negative tail of the range difference is long, although these outliers belong almost exclusively to the case with a 50-m terrain correlation length. A total of 494 differences exceeded 10 m in magnitude. Of these, 469 occurred in the case with a 50-m correlation length, 17 occurred in the case with a 100-m correlation length, and 8 occurred in the case with a 150-m correlation length, evidencing a correlation between range differences and terrain roughness.

In summary, the analytical MRPD range expression in Equation (18) closely matches the empirically found MRPD range values in almost all cases explored in a comprehensive simulation study. Rare exceptions occur more often over rough terrain whose median slope exceeds approximately 10°. Caution is recommended when applying Equation (18) in such cases.

5.4 Jacobian Model Verification

Equations (21) and (26) provide analytical expressions for the position and velocity Jacobians of the MRPD measurement function. Again, these expressions are based on thin Doppler cone and planar terrain assumptions. To verify the validity of these expressions and, by proxy, these two assumptions, the analytical Jacobians were compared against their corresponding empirical finite-difference approximations. The analytical Jacobians were calculated with Equations (21) and (26) using a planar terrain approximation around a once-solved MRPD ground point; the finite-difference Jacobians were approximated by repeatedly perturbing the radar position and velocity and searching for the MRPD ground point. Importantly, because the finite difference was computed over the entire search, this approach takes into account the effect of a non-zero-width Doppler bin and any curvature in the terrain.

The analytical and finite-difference position and velocity Jacobians were compared with two metrics designed to capture angle and magnitude differences between their respective vectors. Let v_α ∈ ℝ³ and v_f ∈ ℝ³ be two candidate analytical and finite-difference Jacobian vectors, respectively, to enable a comparison. The first metric computes the angle Δ θ ∈ [0, +180] in degrees between the two vectors:

$\Delta \theta ={\cos }^{-1}\left[\frac{v_a^T{v}_f}{\Vert v_a\Vert \Vert v_f\Vert }\right]$39

The second metric computes the percent difference in magnitude Δr > 0 between the two vectors:

$\Delta r=\left|\frac{\Vert v_a\Vert -\Vert v_f\Vert }{\Vert v_f\Vert }\right|\times 100$40

In this comparison, the finite-difference Jacobian should be interpreted as the true Jacobian, as the finite difference is applied over the entire search process and takes into account the Doppler bin width and terrain nonlinearities. Of course, the discrete nature of the terrain maps introduces errors in the finite-difference Jacobian; thus, in truth, neither Jacobian is perfectly accurate.

These two metrics were computed for each position and velocity Jacobian associated with each MRPD measurement. Figures 6 and 7 present the CDFs of the position and velocity Jacobian difference metrics, respectively, for each of the five terrain maps. The CDF plots are again split across the five terrain maps to illustrate how the terrain roughness affects the Jacobians. Similar to the range difference CDF plot, there exists a nearly monotonic growth in Jacobian angle and magnitude differences as the terrain correlation length decreases.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Cumulative distributions of the position Jacobian angle (left) and magnitude (right) metrics

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

Cumulative distributions of the velocity Jacobian angle (left) and magnitude (right) metrics

Overall, the analytical and finite-difference Jacobians are similar, although there are some remarkable differences. The analytical position Jacobians closely match the finite-difference position Jacobians in almost all cases: fewer than ≈2% of the position Jacobian vectors differ by more than 10° in direction, and fewer than ≈3% differ by more than 10% in magnitude. The maximum angular difference was 45°, and the maximum magnitude difference was 54%. These differences are remarkably small given the thin Doppler cone and planar terrain approximations invoked in the analytical Jacobian expression.

The analytical velocity Jacobians do not match as closely: while only ≈3% of angles differ by more than 10°, almost 1% differ by more than 170° — a complete sign reversal. Moreover, almost ≈30% of the velocity Jacobians differ by more than 10% magnitude, and some differed by more than 10⁶%. This difference in correspondence between the position and velocity Jacobians can be explained by two factors. First, the magnitude comparison metric in Equation (40) exaggerates the magnitude difference of vectors when they are small, owing to a division by near-zero. The magnitudes of the position Jacobians are bounded below by unity, but no such lower bound exists for the magnitude of the velocity Jacobians, allowing their magnitude differences to become exaggerated in some cases. Second, the ground-projected location of a Doppler bin is sensitive to velocity errors, and significant movement of a Doppler bin can violate planar terrain assumptions when terrain is extremely rough. A perturbation in radar position shifts the location of a Doppler bin by an approximate maximum distance equal to the position perturbation. In contrast, a perturbation in radar velocity changes the Doppler angle, thereby shifting the location of a Doppler bin by an amount proportional to the radar’s velocity and altitude. This angular relationship is captured by the $\partial \alpha /\partial v_r$ expression in Equation (24). For this Monte Carlo simulation, a 10-cm/s velocity perturbation (the velocity finite-difference step size) can shift the ground-projected location of a Doppler bin by ≈17 m — a 34-m span if both perturbation directions are considered. This represents a significant fraction of the terrain correlation length for the 50- and 100-m correlation length maps and likely violates local planar terrain assumptions.

In summary, the analytical MRPD position Jacobians computed from Equation (21) closely match the finite-difference-approximated Jacobians in almost all scenarios. The analytical MRPD velocity Jacobians computed from Equation (26), however, do not match as closely to their finite-difference counterparts and can incur significant angle and magnitude errors when terrain is extremely rough. Caution is recommended when applying Equation (26) in cases where the median terrain slope exceeds approximately 10°.

5.5 Position CRLB Model Verification

Equations (31) and (33)–(35) provide analytical expressions for the FIM of the MRPD measurement function. These expressions yield some insights into the behavior of the TITAN algorithm’s CRLB, yet other properties are hidden. For example, the position block FIM expression in Equation (33) indicates that the primary positioning performance of the algorithm is related to the outer product of the terrain normals at the MRPD ground points and, thus, the variance of these normals. One might assume, then, that the statistics of the terrain normals at the MRPD ground points are equal to the statistics of the terrain normals of the map and that the map statistics could be substituted into Equation (33). However, this is not the case. Indeed, it will be shown that the MRPD process shapes the statistics of the terrain normals at the MRPD ground points. This behavior is not obvious from Equations (31) and (33)–(35) and is explored via Monte Carlo simulations.

The CRLB bounds the expected error covariance of any unbiased estimator. A metric that is closely related to the error covariance is the dilution of precision (DOP). The DOP captures the performance of an estimator along prescribed directions independent of the system measurement noise. Let ${\overline{J}}_N\in {ℝ}^{6\times 6}$ be the FIM given by Equation (38) associated with N scalar MRPD measurements captured at a single measurement epoch. Let us assume that the measurement noise is identical and independent (R = σ²I_N×N). The lower bound on the error covariance P_N is as follows:

$\begin{array}{c}{P}_N:={\overline{J}}_N^{-1} ={\sigma }^{2}E{\left[H{\left[{\hat{O}}^{\prime }\right]}^{T}H\left[{\hat{O}}^{\prime }\right]\right]}^{-1}\\ ={\sigma }^2\left[\begin{array}{cc}{C}_f^{\ell }& 0_{3\times 3}\\ 0_{3\times 3}& C_f^{\ell }\end{array}\right]\left[\begin{array}{ccccc}{\sigma }_{\text{AT}}^2& ·& ·& ·& ·\\ ·& {\sigma }_{\text{CT}}^2& ·& ·& ·\\ ·& ·& {\sigma }_{\text{VT}}^2& ·& ·\\ ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·\\ ·& ·& ·& ·& ·\end{array}\right]{\left[\begin{array}{cc}{C}_f^{\ell }& 0_{3\times 3}\\ 0_{3\times 3}& C_f^{\ell }\end{array}\right]}^T\end{array}$41

where $H\left[{\hat{O}}^{\prime }\right]$ is a random realization of Equation (30). $C_f^{\ell }\in SO\left(3\right)$ is an orthonormal frame transformation matrix that transforms a vector expressed in the flight frame f (along-track, cross-track, vertical) to a vector expressed in the local navigation frame in which the FIM is calculated ℓ (e.g., north, east, down). · are placeholders for terms unimportant to the current discussion, and σ_AT, σ_CT, σ_VT > 0 are the along-track, cross-track, and vertical DOP values, respectively. These DOP values may be interpreted as the expected along-track, cross-track, and vertical positioning error per unit of range measurement noise.

The expectation in Equation (41) was approximated with a simple average over the 500 Monte Carlo iterations for each terrain map. Figure 8 presents the resulting along-track, cross-track, and vertical DOP values versus the number of Doppler bins for each of the five terrain maps. In each case, the DOP values were extracted from P_N, which was calculated in two ways: (1) by inverting the analytical FIM produced via Equations (31) and (33)–(35) and (2) by inverting an evaluation of Equation (31) using Jacobians approximated by a finite-difference approach. Both methods are included to compare the analytical expression with all of its implicit assumptions against an empirical value that accounts for nonlinearities in terrain and non-zero Doppler bin widths. Over all five terrain maps and over all Doppler bins, the analytical and finite-difference-approximated DOP values show strong agreement. This strong agreement indicates that the differences in Jacobian angles and magnitudes shown in the previous section amount to little average difference in fundamental navigation performance.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

Cross-track, along-track, and vertical DOP as a function of the number of Doppler bins across the five terrain maps

Before continuing, readers should take note of the assumptions implicit in and potential limitations of Figure 8. The DOP values depicted in this figure were calculated by averaging instance-specific FIMs for an aircraft flying at an altitude of 9000 ft at a velocity of 80 m/s with a nadir-pointed radar and evaluated over synthetic isotropic terrain. A change in any of these simulation parameters would affect the α, θ, $\hat{o}$, p_r, and v_r variables, which control the resulting DOP values. Furthermore, while prior work has concluded that the performance of the TERCOM algorithm with synthetic terrain closely matches the performance of the algorithm with real terrain when the terrain statistics are similar (Cannon, 1974; Cannon Jr, 1978), such a claim has not been verified for the TITAN algorithm. Verification of this claim is deferred to future work, as it exceeds the scope of this publication. Consequently, readers are cautioned against direct interpretation of the DOP results in Figure 8 and are instead directed to consider the trends in the DOP.

The DOP trends in Figure 8 reveal several remarkable properties of the TITAN algorithm’s CRLB that are perhaps not obvious from Equations (31) and (33)–(35). First, the vertical DOP is always smaller than the along-track DOP, which is always smaller than the cross-track DOP. It is not surprising that the vertical DOP is the smallest, as terrain normals are usually distributed around the vertical. What is surprising, however, is that the along-track DOP is smaller than the cross-track DOP. The FIM expressions indicate no such preference in direction. In fact, the synthetic terrain guarantees that the population statistics of the terrain normals have no preference for direction — Figure 4 indicates that the terrain normals are approximately circularly distributed about the vertical. One may conclude that the MRPD measurement process shapes the statistics of the terrain normals at the MRPD ground points such that there is a greater angular variance in the terrain normals at the MRPD ground points in the along-track direction than in the cross-track direction. Second, whereas the along-track DOP values decrease with increasing terrain roughness, the cross-track DOP values remain nearly stationary across the five terrain maps. This result is also surprising, as one would expect that increasing terrain roughness would increase the variance of the terrain normals in the cross-track direction at the MRPD ground points, thus decreasing the cross-track DOP. Instead, Figure 8 indicates that this cross-track variance is independent of the population variance of the terrain normals.

Both of these remarkable properties can be explained by the same process. Recall that the MRPD measurement function is a constrained optimization problem. In fact, the analytical expression models the MRPD measurement function as an equality-constrained optimization problem. A necessary condition of equality-constrained optimization problems is that the contours of the optimization function and constraints must be tangent at the minimizer — a consequence revealed by Lagrange multipliers. With this property in mind, suppose that the MRPD optimization problem were broken down into two independent optimization problems: (1) a distance minimization problem in the along-track direction subject to terrain and Doppler constraints and (2) a distance minimization problem in the cross-track direction subject only to terrain constraints. This decomposition is motivated by the observation that the Doppler constraint effectively constrains only the location of MRPD ground points in the along-track direction. Of course, the 3D distance minimization problem subject to 3D Doppler and terrain constraints cannot be exactly decomposed into two 2D problems, but let us accept the decomposition as an explanatory mental model, if not a rigorous mathematical model.

Figure 9 illustrates the two decomposed 2D optimization problems. The rear-view case illustrates the cross-track optimization problem wherein the distance minimization problem is subject only to the terrain constraint; the profile-view case illustrates the along-track optimization problem wherein the distance minimization problem is subject to both the terrain and Doppler constraints. In the rear-view case, with only the terrain constraint, the Lagrange multiplier conditions require that the constant-range circle be tangent to the terrain at the minimum-range point. It follows that the cross-track look vector from the radar’s position to the MRPD ground point must be equal and opposite to the cross-track terrain normal vector. Because the MRPD measurement function minimizes distance, the cross-track look vectors (and, therefore, the terrain normal vectors) will tend toward the vertical because distances to the ground tend to increase with look vectors further from the nadir. To summarize, the cross-track terrain normal vectors are locked to the cross-track MRPD look vectors, which are concentrated near the nadir.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

Illustration of the model decomposition of the 3D optimization problem into two 2D optimization problems. The cross-track optimization problem is only subject to terrain constraints, whereas the along-track optimization problem is subject to both terrain and Doppler constraints. This difference in constraints results in different terrain normal vectors at the MRPD ground points.

In the profile-view case, both the terrain and Doppler constraints are present. That is, in the along-track direction, the 2D optimization problem is subject to two constraints. In instances where the number of constraints is equal to the dimensionality of the problem (and provided a unique solution exists), no optimization is possible, as the constraints alone determine the location of the minimizer. Accordingly, no constraints are placed on the Jacobians of the constraint functions (terrain normal vectors). For the MRPD optimization problem, this means that the along-track terrain normal vectors at the MRPD ground points can be arbitrarily oriented. Figure 9 depicts the arbitrary nature of these terrain normals — some terrain normals at MRPD ground points are aligned with the MRPD look vector whereas others are nearly orthogonal. To summarize, no constraints are placed on the along-track terrain normal vectors; thus, they are allowed to vary with the population statistics of the terrain.

Supposing that the optimization problem decomposition is reasonable, these two observations — (1) an enforced relationship between cross-track look vectors and cross-track terrain normals and (2) no such relationship in the along-track direction — explain the two remarkable properties of DOP trends shown in Figure 8. First, in instances where the terrain normals vary more than the look vectors, the cross-track DOP is expected to be larger than the along-track DOP, as the variance of the terrain normals in the cross-track direction is attenuated down to the look vector variance, whereas no such attenuation occurs for the along-track normals. This difference in terrain normal variance explains the difference in observed along-track and cross-track DOP values within each terrain map case. Second, because all five synthetic terrain maps were constructed with the same standard deviation for terrain elevation (σ_h), the synthesized terrain is, on average, no closer to the radar in any particular case. As a result, the horizontal distribution of the nearest points in the cross-track direction — and, thus, the cross-track look vectors — is similar across all cases. Because the cross-track terrain normals are locked to the cross-track look vectors, this result implies that the distribution of the cross-track terrain normals is also invariant across the five terrain maps, explaining the observed independence between the cross-track DOP and terrain roughness.

These conclusions are, of course, conditioned on the validity of the 2D decomposition of the 3D optimization problem, which must be evidenced. Let $\hat{o}={\left[\begin{array}{ccc}{x}_{\hat{o}}& y_{\hat{o}}& z_{\hat{o}}\end{array}\right]}^T$ be the terrain normal vector at the MRPD ground point, and let $\hat{r}={\left[\begin{array}{ccc}{x}_p& y_p& z_p\end{array}\right]}^T$ be the look vector from the radar to the MRPD ground point. Suppose that the radar is flying in the x direction, that the terrain normal is pointed up toward the radar, and that the look vector is pointed down toward the terrain. We define the along-track and cross-track angles of the terrain normals as follows:

${\Psi }_{\text{AT}}={\tan }^{-1}\left[\frac{x_{\hat{o}}}{z_{\hat{o}}}\right]$42

${\Psi }_{\text{CT}}={\tan }^{-1}\left[\frac{y_{\hat{o}}}{z_{\hat{o}}}\right]$43

Similarly, we define the along-track and cross-track angles of the MRPD look vectors as follows:

${\Phi }_{\text{AT}}={\tan }^{-1}\left[-\frac{x_{\hat{r}}}{z_{\hat{r}}}\right]$44

${\Phi }_{\text{CT}}={\tan }^{-1}\left[-\frac{y_{\hat{r}}}{z_{\hat{r}}}\right]$45

If the 2D decomposition model is reasonable, then it follows that the cross-track angles of the terrain normals and look vectors will be nearly equal, whereas the along-track angles of the terrain normals and look vectors will be uncorrelated. Figure 10 displays the distribution of the difference in angles between the terrain normals and look vectors over all of the Monte Carlo simulations. The cross-track angles are strongly concentrated — fewer than 5% of the angles differ by more than 1°. In contrast, the along-track angles have a wide distribution — more than 90% of the angles differ by more than 1°. This alignment of cross-track angles and lack of alignment of the along-track angles support the conclusions drawn from the modeled 2D decomposition of the optimization problem and indicate that the model is reasonable.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 10

Distribution of the difference in marginal along-track and cross-track angles between the terrain normals at the MRPD ground points and the look vectors to the MRPD ground points. The small difference in angles in the cross-track direction indicates strong alignment between the normal vectors and look vectors in the cross-track direction.

The findings from the CRLB verification can be summarized as follows:

1. There is excellent agreement between the empirical finite-difference position CRLB and the position CRLB derived from the FIM expressions in Equations (31) and (38).

2. The vertical DOP is smaller than the along-track DOP, which is smaller than the cross-track DOP.

3. The difference in cross-track and along-track DOP values is related to the fact that the Doppler constraint effectively acts only in the along-track direction.

4. The cross-track DOP is related to the variance of the radar look vectors, whereas the along-track DOP is related to the variance of the terrain normals.