Mapping Target Location from Doppler Data

TL;DR

This paper introduces an algorithm that determines a target location curve on Earth's surface using Doppler shift measurements from UAVs or satellites, incorporating terrain data and environmental factors.

Contribution

The novel algorithm combines Doppler data, terrain mapping, and environmental considerations to accurately locate a stationary emitter on Earth's surface.

Findings

Successfully maps target location curves using Doppler data.

Incorporates Earth's rotation and atmospheric effects into the model.

Provides a method to handle measurement errors and uncertainties.

Abstract

In this paper, we present an algorithm for determining a curve on the earth's terrain on which a stationary emitter must lie according to a single Doppler shift measured on an unmanned aerial vehicle (UAV) or a low earth orbit satellite (LEOS). The mobile vehicle measures the Doppler shift and uses it to build equations for a particular right circular cone according to the Doppler shift and the vehicle's velocity, then determines a curve consisting of points which represents the intersections of the cone with an ellipsoid that approximately describes the earth's surface. The intersection points of the cone with the ellipsoid are mapped into a digital terrain data set, namely Digital Terrain Elevation Data (DTED), to generate the intersection points on the earth's terrain. The work includes consideration of the possibility that the rotation of the earth could affect the Doppler shift, and of the errors resulting from the non-constant refractive index of the atmosphere and from lack of precise knowledge of the transmitter frequency.

Equations

$$\begin{array}{rl}x& =\left(\frac{a}{\chi }+h\right)\cos \varphi \cos \lambda \\ y& =\left(\frac{a}{\chi }+h\right)\cos \varphi \sin \lambda \\ z& =\left(\frac{a(1-e^2)}{\chi }+h\right)\sin \varphi \end{array}$$

$$\chi =\sqrt{1-e^2{\sin }^2\varphi }$$

$$\begin{array}{rl}& \lambda =\arctan \left(\frac{y}{x}\right)\\ & \varphi =\arctan \left(\frac{z(1-f)+(2f-f^2)a{\sin }^3\mu }{(1-f)(p-(2f-f^2)a{\cos }^3\mu )}\right)\\ & h=p\cos \varphi +z\sin \varphi -a\sqrt{1-(2f-f^2)si{n}^2\varphi }\end{array}$$

$$\mu =\arctan \left(\frac{z}{p}(1-f)+\frac{(2f-f^2)az}{rp}\right)$$

$$\displaystyle\left[\begin{array}[]{c}de\\ dn\\ du\end{array}\right]$$

$$=$$

$${\mathbf{R}}_{\text{Body}}^{\text{ENU}}=\left[\begin{array}{ccc}\sin \gamma \cos \beta & \cos \alpha \cos \gamma +\sin \alpha \sin \gamma \sin \beta & -\sin \alpha \cos \gamma +\cos \alpha \sin \gamma \sin \beta \\ \cos \gamma \cos \beta & -\cos \alpha \sin \gamma +\sin \alpha \cos \gamma \sin \beta & \sin \alpha \sin \gamma +\cos \alpha \cos \gamma \sin \beta \\ \sin \gamma & -\sin \alpha \cos \beta & -\cos \alpha \cos \beta \end{array}\right]$$

$$h(\varphi ,\lambda )=H(\varphi ,\lambda )+N(\varphi ,\lambda )$$

$$\cos \psi =(\delta /f_0)/(|\mathbf{v}|/c)$$

$$f_r=f_0\left(1-\frac{1}{c}\left[\frac{d(\mathbf{p}-\mathbf{r})}{dt}\right]\cdot \frac{\mathbf{p}-\mathbf{r}}{\Vert \mathbf{p}-\mathbf{r}\Vert }\right)$$

$$\frac{d{\mathbf{e}}_i}{dt}=\mathbf{\Omega }\times {\mathbf{e}}_{i}i\in \{x,y,z\},$$

$$\frac{d\left(\mathbf{p}-\mathbf{r}\right)}{dt}$$

$$f=f_0\left(1-\frac{1}{c}\left(\sum _i\frac{d\left(p^i-r^i\right)}{dt}{\mathbf{e}}_i\cdot \frac{\mathbf{p}-\mathbf{r}}{\Vert \mathbf{p}-\mathbf{r}\Vert }\right)\right)$$

$$\frac{x^2}{d^2}+\frac{y^2}{d^2}-z^2=0$$

$$\frac{(x-r_x{)}^2}{d^2}+\frac{(y-r_y{)}^2}{d^2}-(z-r_z{)}^2=0$$

$$\mathbf{R}=\left[\begin{array}[]{ccc}\frac{\alpha\gamma}{(\alpha^{2}+\beta^{2})^{1/2}}&-\frac{\beta}{(\alpha^{2}+\beta^{2})^{1/2}}&\alpha\\ \frac{\beta\gamma}{(\alpha^{2}+\beta^{2})^{1/2}}&\frac{\alpha}{(\alpha^{2}+\beta^{2})^{1/2}}&\beta\\ -(\alpha^{2}+\beta^{2})^{1/2}&0&\gamma\end{array}\right]$$

$$\mathbf{R}(\theta,\gamma)=\left[\begin{array}[]{ccc}\gamma\cos\theta&-\sin\theta&\sqrt{1-\gamma^{2}}\cos\theta\\ \gamma\sin\theta&\cos\theta&\sqrt{1-\gamma^{2}}\sin\theta\\ -\sqrt{1-\gamma^{2}}&0&\gamma\end{array}\right]$$

$$\mathbf{R}(\theta,1)=\left[\begin{array}[]{ccc}\cos\theta&-\sin\theta&0\\ \sin\theta&\cos\theta&0\\ 0&0&1\end{array}\right]$$

$$\frac{{\bar{x}}^2}{d^2}+\frac{{\bar{y}}^2}{d^2}-{\bar{z}}^2=0$$

$${\bar{\mathbf{p}}}^{\top }\bar{\Lambda }\bar{\mathbf{p}}=0$$

$$\bar{\Lambda }=\text{mbox}diag[d^{-2},d^{-2},-1]$$

$$\left[\begin{array}[]{c}x\\ y\\ z\end{array}\right]=\mathbf{R}\left[\begin{array}[]{c}\bar{x}\\ \bar{y}\\ \bar{z}\end{array}\right]$$

$${\mathbf{p}}^{\top }\mathbf{R}\bar{\Lambda }{\mathbf{R}}^{\top }\mathbf{p}=0$$

$$\left[x-r_{x}\;y-r_{y}\;z-r_{z}\right]\mathbf{R}\operatorname{diag}[d^{-2},d^{-2},-1]\mathbf{R}^{\top}\left[\begin{array}[]{c}x-r_{x}\\ y-r_{y}\\ z-r_{z}\end{array}\right]=0$$

$$\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1$$

$$\begin{array}{r}x=r_x+{\alpha }_{c}s\\ y=r_y+{\beta }_{c}s\\ z=r_z+{\gamma }_{c}s\end{array}$$

$$\begin{array}{rl}{\alpha }_r& =\cos (\zeta )\cos (\eta )\\ {\beta }_r& =\cos (\zeta )\sin (\eta )\\ {\gamma }_r& =\sin (\zeta )\end{array}$$

$${\mathbf{d}}_t=\mathbf{R}{\mathbf{d}}_r.$$

$$a_s{s}^2+b_{s}s+c_s=0$$

$$\begin{array}{rl}{a}_s& =q_1{\alpha }_c^2+q_2{\beta }_c^2+q_3{\gamma }_c^2\\ b_s& =2q_1{r}_x{\alpha }_c+2q_2{r}_y{\beta }_c+2q_3{r}_z{\gamma }_c\\ c_s& =q_1{r}_x^2+q_2{r}_y^2+q_3{r}_z^2+q_0;\end{array}$$

$$D_s=b_s^2-4a_s{c}_s$$

Full text

I Introduction

This paper is concerned with the following scenario. A stationary emitter is located at an unknown location on the surface of the earth. The frequency of its transmissions is known. An unmanned aerial vehicle (UAV) or a low earth orbit satellite (LEOS) receives the transmission. Using the Doppler shift, it is required to determine a curve, depicted using a 3D map tool, on which the transmitter must lie.

An idealized version of this scenario might assume that the earth is a sphere, or possibly an ellipsoid, that all measurements (of position, velocity and frequency) were noise free and that the rotation of the earth can be neglected. In this case, the UAV or LEO receiver data serves to determine a right circular semi-infinite cone on which the emitter must lie, the axis coinciding with the velocity vector direction, the semi-angle derived from the Doppler shift, and the apex obtained from the receiver position. The intersection of this cone with the sphere/ellipsoid corresponding to the earth then provides the solution to the problem.

There are however a number of complicating factors. These include

the multiple related but differing representations of the earth surface, including the WGS84 ellipsoid, defined by a constant gravitational equipotential surface; the EGM96 geoid, as an approximate ellipsoid with fine structure topographical variations capturing the height above sea level of points on the earth’s surface [1]. 2. 2.

the different types of intersection/non-intersection that can occur between a cone and a ellipsoid. There may be no intersection; there may be a tangency; there may be a single curve; there may be two separate, i.e. non-intersecting, curves (one of which may be in the ’shadow‘ of the earth as seen from the receiver), etc. 3. 3.

the data sources available to represent the earth, e.g. the parameters of the WGS84 model [1], the DTED data base [2], which indicates for samples spaced at a known interval the height of a point on the surface of the earth above the geoid. DTED is a matrix of terrain elevation values which provides basic quantitative data for systems and applications that require terrain elevation, slope, and/or surface roughness information. DTED data is uniformly spaced in angle (not distance). 4. 4.

the fact that the earth is rotating about a polar axis, which means that given two points that are stationary with respect to the earth’s surface except for the two poles, there is actually a relative velocity; roughly speaking this is a Coriolis effect. The value depends on the position of the two points. In consequence, it would seem that the determination of a cone using a Doppler shift, which amounts to indirect use of velocity data of the receiver (and assuming a stationary emitter), ought to take into account this relative velocity. 5. 5.

the fact that a nominal rather than exact value for the emitter frequency may be known. (By way of example, a mobile phone tower may will an precisely known value, but a mobile phone may only have a nominal value). When a nominal value is used, there will be an error in the computation of Doppler shift and therefore of the cone and it is intersection with the Earth’s surface is moved. 6. 6.

the fact that transmissions between a ground-located transmitter and a LEO receiver pass through a sufficiently long distance that an assumption of uniform value of the refractive index, or equivalently straight line propagation, is dangerous. (In contrast, if a UAV rather than LEO receives the transmission, the assumption of straight line propagation is almost certainly reasonable.) The refractive index varies with altitude. This variation will be significant between a ground receiver and a satellite, but minor between a ground receiver and a low-altitude UAV.

A significant part of the contents of this paper seeks to deal, at least to some degree, with these issues.

Our work is of course not the first work on this subject. Prior contributions, some with more idealized assumptions than those we seek to have, certainly exist. In [3, 4, 5], geolocation of radio frequency (RF) emitters using Doppler shift measurements from mobile vehicles was considered. These works investigated scenarios assuming the vehicle path and the emitter are located in a plane, and so without considering the shape of the earth’s surface. In [6], a set of solutions from the combination of time difference of arrival (TDOA) and frequency difference of arrival (FDOA) measurements for localizing an emitter with known altitude above the earth surface was proposed. However, the method assumed an ellipsoid earth model without considering any height information. This restriction also appeared in [7, 8], which had the similar problem settings to those in [6].

The structure of the paper is as follows. In the next section we provide a general review of various coordinate systems used for representing points on or above earth’s surfaces, in this paper, relating them to the World Geodetic System (WGS) 1984 , and its related Geoid system for representing the earth’s terrain. We introduce Digital Terrain Elevation Data (DTED), which is used as a data reference to index specific points on the earth’s terrain. We also review standard coordinate systems that represent points on or above the earth surface, including earth-centered-earth-fixed coordinate systems, geographic coordinate systems and body coordinate systems, and explain their connections. In section III, we present fundamental Doppler-shift equations and identify two possible sources of Doppler shift in our scenario, viz. the motion relative to the earth’s surface of the vehicle and the effect of the earth’s rotation, and reveal the fact that the earth’s rotation actually does not affect the Doppler shift, due to the fulfillment of a certain orthogonality condition. Though in our scenario, we consider a stationary emitter, the conclusion that the earth’s rotation does not affect the Doppler shift is not dependent on this staitonarity, as the later derivation shows. We present the equations to build a particular right circular cone with the knowledge of the sensing vehicle position, velocity and the measured Doppler shift. In Section IV-A, we present an algorithm to find the intersections of a right circular cone with the WGS84 ellipsoid and identify different types of intersections. In Section IV-B, we present our developed method to find the cone-earth-terrain intersection curve where the emitter lies, by using the formerly built cone-earth-ellipsoid intersections and the DTED. We demonstrate the effectiveness of this method with some examples 111Before moving to the next section, we refer the readers to the appendix for the notation and abbreviations used throughout this paper..

II Preliminaries

II-A WGS84 ellipsoid and EGM96 geoid

In general, global geodetic applications require three different surfaces to be clearly defined:

The earth’s topographic surface, which includes the landmass topography and the ocean bottom topography 2. 2.

A geometrical or mathematical reference surface, which is an ellipsoid with known semi-axes. 3. 3.

An equipotential surface, called the geoid. (Potential refers to gravity). Because of height variation, and because of density variation within the solid earth, this is bumpy rather than a smooth ellipsoid, and it does not coincide with the topographic surface.

The World Geodetic System (WGS) is a standard for use in cartography, geodesy, and satellite navigation including GPS. The latest revision is WGS84(G1762), established in 1984 and last revised in 2013. It comprises a standard coordinate system for the Earth, a standard ellipsoidal reference surface (the datum or reference ellipsoid) for raw altitude data, and a gravitational equipotential surface (the geoid) that defines the nominal sea level.

The coordinate system defined in WGS84 is depicted in the next subsection.

The WGS84 reference ellipsoid is a mathematically defined surface that approximates the truer shape of the Earth. It is used as a preferred surface on which geodetic computations are performed and point coordinates such as latitude, longitude, and elevation are defined. The ellipsoid is defined in terms of the semi-major axis value $a=6378137.0\text{m}$ and a flattening coefficient $f=1/298.257223563$ (from which the semi-minor axis value $b=6356752.314245\text{m}$ can be determined). The fact that an ellipsoid rather than a sphere is involved gives rise to an important observation. Whereas in normal life, the notion of a vertical direction is considered to be identical with the notion of the direction of the force exerted by gravity on a body on or above the earth’s surface (disregarding sign), these two directions cannot strictly be identified.

The geoid is a smooth but highly irregular surface and has been built as a mathematical representation of the surface of the earth’s gravity field. The geoid is widely used to describe mean sea level (MSL) since the surface of the gravity field coincides approximately with the mean sea level. Disregarding the terrain elevation on continents, the geoid is a much more accurate description of the true physical shape of the earth, than the WGS84 ellipsoid, though it is not of course identical. The WGS84 reference ellipsoid is the baseline of EGM96 geoid, i.e., the geoid undulations are with respect to the WGS84 ellipsoid. EGM96 geoid is a data grid overlay onto the WGS84 ellipsoid and can be used to find MSL for almost any point on the Earth at any given latitude and longitude.

II-B Coordinate systems

In this subsection, we review standard material concerning the different coordinate systems used for representing points on or above the earth’s surface

The first of these systems is an earth-centered, earth-fixed (ECEF) cartesian coordinate system. A stationary point on the surface of the earth has constant ECEF coordinates. The origin is taken to be the earth’s center of mass. The $z$-axis passes from the origin through a point on the surface termed the reference pole (effectively the north pole), and the $x$-axis passes through a zero meridian (effectively the longitudinal line through Greenwich, UK). The $y$-axis is chosen to ensure a right-handed orthogonal coordinate system.

The second coordinate system is the geographic coordinate system, which enables every location on Earth to be specified by a triple consisting of latitude $\phi$, longitude $\lambda$ and height $h$. The latitude of a point on or above the earth’s surface (ellipsoid surface) is the angle between the equatorial plane and the normal at that point. Note that the normal at a point on an ellipsoid surface does not pass through the centre, except for points on the equator or at the poles. The height of a point is the vertical distance from the point to some surfaces (e.g. WGS84 ellipsoid or EGM96 geoid).

The third coordinate system is the body coordinate system, which is directly defined on the body of the vehicle. Its origin is located at the center of gravity of the vehicle. The $x$-axis of the body coordinate system points froward; the y-axis points to the right of the $x$-axis, perpendicular to the $x$-axis; the $z$-axis points down through the bottom the vehicle, perpendicular to the $x$-$y$ plane.

The last coordinate system is the vehicle-carried navigation coordinate system, which is usually viewed as recording east, north and up displacements from some agreed point of origin. We denote the vehicle-carried navigation coordinate system with the term ENU coordinate system in this paper. Its origin is defined as the center of the gravity of the vehicle as well. A tangent plane is fitted to this fixed origin point, with the east, north and up axes pointing in the obvious directions.

II-C Connection between coordinate systems

Navigation instruments such as the global positioning system (GPS) and the inertial navigation system (INS) may well determine location (latitude, longitude and height) and attitude (roll, pitch, yaw) in relation to some agreed vehicle and some system is needed to marry the use of measurements derived from such systems and ECEF measurements.

We first consider the connection between ECEF and geographic coordinate systems. The output of a GPS receiver is latitude $\phi$, longitude $\lambda$ and height $h$ in the geographic coordinate system, denoted by $(\phi,\lambda,h)$ (the vertical baseline is the WGS84 ellipsoid) [9]. The relevant equations for converting $(\phi,\lambda,\text{h})$ to $(x,y,z)$ in the ECEF coordinate system are [10]:

[TABLE]

where $e^{2}=6.69437999014\times 10^{-3}$ is the square of the first eccentricity, and

[TABLE]

The inverse transformation from ECEF coordinates to geographic coordinates can be described by the following equations [11]:

[TABLE]

where $p=\sqrt{x^{2}+y^{2}}$, $r=\sqrt{p^{2}+z^{2}}$, $f$ is the flattening of the ellipsoid and $\mu$ is a parameter calculated according to

[TABLE]

We now present the connections between ECEF and ENU coordinate systems. Suppose that the origin for ENU coordinates is defined by a point $(x,y,z)$ in ECEF coordinates. Suppose that a point near to the specified point is defined by the coordinate differences $(dx,dy,dz)$ between the coordinates of the nearby point, viz, $(x+dx,y+dy,z+dz)$ and the coordinates of the specified origin point $(x,y,z)$. The ENU coordinates of the nearby point are $(de,dn,du)$, the origin point of course having ENU coordinates $(0,0,0)$. The question arises: how is $(de,dn,du)$ related to $(dx,dy,dz)$. This can be answered by understanding that the orientation of ENU coordinates is determined by rotating the ECEF coordinates; the first rotation is about the $z$-axis, by $\lambda$ degrees (corresponding to the longitude), and then rotating about the new $x$-axis (obtained from the old $x$-aixs through rotation of $\lambda$ degrees in the $x$-$y$ or equatorial plane) by $\phi$ degrees. Consequently we have

\displaystyle\left[\begin{array}[]{c}de\\ dn\\ du\end{array}\right] $\displaystyle=$ \displaystyle\left[\begin{array}[]{ccc}1&0&0\\ 0&-\sin\phi&\cos\phi\\ 0&\cos\phi&\sin\phi\end{array}\right]\left[\begin{array}[]{ccc}-\sin\lambda&\cos\lambda&0\\ \cos\lambda&\sin\lambda&0\\ 0&0&1\end{array}\right]\left[\begin{array}[]{c}dx\\ dy\\ dz\end{array}\right]

(16)

$\displaystyle=$ \displaystyle\left[\begin{array}[]{ccc}-\sin\lambda&\cos\lambda&0\\ -\sin\phi\cos\lambda&-\sin\phi\sin\lambda&\cos\phi\\ \cos\phi\cos\lambda&\cos\phi\sin\lambda&\sin\phi\end{array}\right]\left[\begin{array}[]{c}dx\\ dy\\ dz\end{array}\right]

(23)

We now state the connections between the vehicle-carried navigation system and the body coordinate system. The output of INS is the Euler angles, the set of roll $\alpha$, pitch $\beta$, raw $\gamma$, with respect to the vehicle-carried navigation coordinate system. The rotation matrix $\mathbf{R}^{\text{ENU}}_{\text{Body}}$ from the ENU coordinate system to the body coordinate system can be described by using the Euler angles, with the formula described by (5) at the top of the next page [12]:

II-D Dealing with height data

Height data is available publicly and has been specified precisely in what is known as DTED (Digital Terrain Elevation Data) [13]. As noted earlier, DTED is a digital elevation model, consisting of a matrix of terrain elevation values, corresponding to the height of the ground above the geoid. The DTED format has three levels, termed level 0, level 1 and level 2, corresponding respectively to spacing of approximately 900 meters, 90 meters and 30 meters. These line spacings correspond to specific arcsecond changes in latitude and longitude.

DTED data is evidently indexed to specific points on the earth’s terrain. The horizontal datum is referenced to the WGS84 ellipsoid and the vertical datum is referenced to the EGM96 geoid. For a particular point above the ground, there are three types of height: orthometric height is defined as the height of the terrain above the geoid: ellipsoid height is defined as the height of the terrain above the WGS84 ellipsoid; geoid height is defined as the height of the geoid above the WGS84 ellipsoid. The relationship of orthometric height, ellipsoid height and geoid height is depicted in Fig. 2, where an up arrow indicates a positive value and down arrow indicates a negative value. It is straightforward to see that the orthometric height at a given point with latitude $\phi$ and longitude $\lambda$ can be expressed by

[TABLE]

II-E Problem statement

Our task is one of passing from a set of measurements to a presentation in a useful (visual) form of information provided by those measurements. More specifically, our task is to define on a representation such as Google map a curve corresponding to those points consistent with a single FDOA measurement taken by an aircraft.

This task can be broken down as follows:

Establish the equation of a right circular cone with fixed apex (corresponding to vehicle position), axis (corresponding to vehicle heading) and semi-angle (corresponding to FDOA measurement) in the ECEF coordinate system. 2. 2.

Establish a procedure for determining an array from which the curve of intersection of the two surfaces (the FDOA cone established in the first step and the WGS84 ellipsoid) can be constructed. 3. 3.

Relate the ECEF coordinates to geographic coordinates. 4. 4.

Make an adjustment to cope with ground height above sea level according to DTED.

III Doppler shift and earth rotation

III-A Locus of points of constant Doppler

We now consider the problem of defining the possible emitter locations corresponding to a particular measured Doppler shift. The scenario is that there is a transmitter on the surface of the earth, and a receiver above the earth that can receive signals.

Throughout this section, we make several explicit assumptions. However, with warning, we will relax specific assumptions at certain points. The assumptions are:

The earth’s surface is modelled as the WGS84 ellipsoid 2. 2.

A transmitter is located on the surface of the earth and a receiver mounted in a sensing vehicle is located above the surface of the earth, with the transmitter being visible from the receiver. 3. 3.

The transmitter is stationary, while the sensing vehicle is moving. 4. 4.

The velocity of light is constant and known; as a consequence, straight line propagation occurs. (This assumption is in a sense one of the least justifiable, and its relaxation will be explored later.) 5. 5.

The exact position and velocity of the sensing vehicle are known at the time it measures a Doppler-shifted signal. 6. 6.

The sensing vehicle can infer the Doppler shift associated with a received signal, because it has precise knowledge of the unshifted frequency at the transmitter.

It is straightforward to make a general statement about the nature of the locus of the points corresponding to a particular Doppler shift. Knowledge of the Doppler shift associated with a received signal implies knowledge of a circular cone on which the transmitter must lie. The cone’s apex is at the sensing vehicle position, and axis aligned with the velocity vector of the sensing vehicle. The direction of the axis can be either the same or reverse of the velocity direction, depending on the sign of the measured Doppler shift. If the measured Doppler shift is positive, then the axis of the cone is defined with the same direction of the velocity vector, otherwise the direction of axis is the reverse of the velocity vector, see Fig. 3. The semi-angle, call it $\psi$, can be computed from knowledge of the magnitude of the velocity $|\mathbf{v}|$ of the sensing vehicle, the speed of light $c$, the unshifted transmitter frequency $f_{0}$, and the Doppler shift $\delta$: 222The transmitter is assumed to be stationary, and for the moment we neglect any contribution to the relative velocity of the receiver with respect to the transmitter due to the earth’s rotation. This will be dealt with subsequently.

[TABLE]

The fractional Doppler shift equals the fraction relative to the speed of light of the component of the sensing vehicle’s velocity along the direction to the transmitter.

The FDOA cone and WGS84 ellipsoid are both two-dimensional surfaces in a 3D space. They intersect in at least one point, namely the point where the transmitter is located. Generically, if there is a nonempty intersection, it will be a one-dimensional smooth set, a curve in space in fact. Exceptionally, as a kind of limiting case, the intersection may be a single point. This would arise if the tangent plane to the ellipsoid coincided with the tangent plane to the cone at the common point of intersection.

Our ultimate goal is to the intersection curve.

III-B Coriolis Effect

In this subsection, we highlight the need to consider the fact that the earth is rotating, when considering the determination of the cone associated with a Doppler shift. We demonstrate that the operative kinematics constraints imply that the relative velocity component between an emitter and a receiver due to the earth’s rotation makes zero contribution to the Doppler shift.

Viewed from a point in a coordinate frame that is not earth fixed, but which views the earth as rotating, any two points on the earth’s surface (other than the poles) can be seen to have velocities associated with the rotation of the earth about its axis. Further, any two distinct points excluding the two poles have different velocities on this account. (The yearly motion of the earth around the sun is discounted in making these observations.) It follows that any two points on (or above) the earth’s surface which in an ECEF coordinate frame is concerned would be considered stationary actually have a relative velocity (whose magnitude can be substantial, e.g. some tens of m/sec), and thus potentially gives rise to a Doppler shift between the transmit/receive frequency of an RF signal propagating from one point to the other. Further, if one of the points is fixed on or above the earth’s surface, and the other is moving (relative to the earth’s surface), the relative velocity can be regarded as coming from two causes, the motion relative to the earth’s surface of one point, and the effect of the earth’s rotation (‘Coriolis effect’). It is because a coordinate frame fixed to the earth is actually rotating, albeit with a fixed angular velocity, rather than just translating at a uniform velocity, that this Coriolis effect arises.

We now establish how the relative velocity can be determined, and seek to determine the associated Doppler shift. For this purpose, it is the case to attribute the Coriolis component to rotation of an earth-centred earth-fixed coordinate basis.

In ECEF coordinates we denote the emitter coordinates by using a vector $\mathbf{p}:=p^{x}\mathbf{e}_{x}+p^{y}\mathbf{e}_{y}+p^{z}\mathbf{e}_{z}$ and the receiver coordinates by $\mathbf{r}:=r^{x}\mathbf{e}_{x}+r^{y}\mathbf{e}_{y}+r^{z}\mathbf{e}_{z}$. Where $\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}$ denote unit vectors aligned with the ECEF coordinate axes. Due to the rotation of the earth, these vectors are not stationary, when seen from an inertial frame. The transmitter frequency $f_{r}$ measured at the receiver is

[TABLE]

In computing the derivative, we must allow for motion relative to the ECEF system (values of the coordinates of $\mathbf{p},\mathbf{r}$ change) and, separately, motion of the coordinate axes.

In the ECEF coordinate system the basis vectors rotate so that [14]

[TABLE]

where $\Omega$ is the magnitude of the angular velocity. Consequently the relative velocity between the emitter and the receiver is the sum of the coordinate velocity and the Coriolis effect:

[TABLE]

The second term on the right hand side of (28) is due to the Coriolis effect. Note that this vector is perpendicular to $\mathbf{p}-\mathbf{r}$ and hence doesn’t contribute to the Doppler effect:

[TABLE]

This is exactly the formula we would write down if we failed to take into account the rotation of the earth, or put another, we can simply neglect the effect associated with rotation of the earth.

III-C Equation of an arbitrary right circular cone

In this section, our aim is to find an equation representing a right circular cone with apex coordinates $(r_{x},r_{y},r_{z})$ and unit vector in the direction of velocity given by $(\alpha,\beta,\gamma)$, i.e. subject to $\alpha^{2}+\beta^{2}+\gamma^{2}=1$, in ECEF coordinate system. Note that a mobile vehicle can determine its attitude $(\alpha,\beta,\gamma)$ with the help of the INS devices. The velocity direction vector $(\alpha,\beta,\gamma)$ cannot be directly obtained but can be calculated by transforming a unit vector $(1,0,0)$ from the vehicle’s body coordinate system to the ECEF coordinate system by using the inverse of the rotation matrices provided in (23) and (5).

Our starting point is that, as is well-known and indeed easily checked, (the surface of) a cone with axis corresponding to the $z$-axis and apex at the origin is given by an equation of the form

[TABLE]

The semi-angle is $\tan^{-1}d$.

Our immediate goal now is to understand how to handle translation of the apex and rotation of the axis of the cone. If the apex of the cone is at $(r_{x},r_{y},r_{z})$, then the above equation is replaced by

[TABLE]

To understand how to handle rotation, we first observe a simple Lemma.

Lemma 1**.**

Let $\alpha,\beta,\gamma$ be a set of direction cosines of a real 3-vector, i.e. $\alpha^{2}+\beta^{2}+\gamma^{2}=1$. Then the following matrix is an orthogonal rotation matrix:

[TABLE]

By way of an outline proof, we observe first that it is easily verified that each of the columns is a vector of length 1, and the columns are mutually orthogonal. Further, a straightforward calculation shows that the determinant is 1, assuring that the matrix is a rotation matrix.

When $\alpha=\beta=0$, four entries of the matrix are not well-defined, and so the continuity of the matrix comes into question. However, let $\theta$ be such that $\alpha=\sqrt{1-\gamma^{2}}\cos\theta$ and $\beta=\sqrt{1-\gamma^{2}}\sin\theta$. Note that given the direction cosines $\alpha,\beta$ and $\gamma$, such a $\theta$ always exists and is unique. Let us regard the matrix $R$ as a function of $\theta$ and $\gamma$. Then we can write

[TABLE]

Now if $\alpha,\beta$ tend continuously to zero while obeying $\alpha^{2}+\beta^{2}<0$ (except in the limit), i.e. $\gamma<1$ except in the limit, and if they tend to zero in such a way that $\beta/\alpha$ also approaches a limit, it is evident that $R(\theta,\gamma)$ will approach the limit

[TABLE]

Now suppose a right-circular cone with apex at the origin is such that the axis of the cone has direction cosines $\alpha,\beta,\gamma$. Assume temporarily that a coordinate basis with coordinates $\bar{\mathbf{p}}=[\bar{x}~{}\bar{y}~{}\bar{z}]^{\top}$ is established with the same origin and with the $\bar{z}$-axis coinciding with the axis of the cone. The equation of the cone in the new coordinate basis is given by

[TABLE]

for suitably chosen $d$. Write this equation as

[TABLE]

where

[TABLE]

It is clear that we want correspondence between the line $\bar{x}=0,\bar{y}=0$, i.e. the axis of the cone in $\bar{x},\bar{y},\bar{z}$ space and the axis of the cone in $x,y,z$ space, which is defined by direction cosines $\alpha,\beta,\gamma$. Because of the structure of the last column of $\mathbf{R}$, this correspondence is assured if we take $\mathbf{p}=\mathbf{R}\bar{\mathbf{p}}$ or

[TABLE]

Equivalently, the equation in $x,y,z$ space of the cone follows from (35) as

[TABLE]

The matrix $\mathbf{R}(\theta,1)$ is relevant in a special case. If $\alpha=0,\beta=0,\gamma=1$ define the directions in which the axis of the cone should point, then it is evident that the original cone, with $z$-axis coinciding with the direction of the cone is already correctly positioned. Thus one might reasonably suppose that the transforming rotation matrix should be the identity matrix. The formula for $\mathbf{R}(\theta,1)$ defines a rotation matrix which leaves that axis invariant, as it causes rotation only of the $x,y$ coordinate plane about the $z$-axis, which is normal to it.

Combining the results on translation and rotation, we have the following Lemma.

Lemma 2**.**

Consider a right circular cone with apex at $\mathbf{p}=[r_{x},r_{y},r_{z}]^{\top}$, with axis given by the unit norm vector of direction cosines $[\alpha,\beta,\gamma]^{\top}$ and with semi-angle given by $\tan^{-1}d$. Then the equation of this cone is

[TABLE]

Our next major task is to study the intersection of this conical surface with the ellipsoid defined by the earth.

IV Determining the intersection of the FDOA cone with the earth’s terrain

Here is a summary of what is required.

The equation of WGS84 ellipsoid (with the positive $z$-axis identified with a line starting at the center of the earth and passing through the north pole) in ECEF coordinates is

[TABLE]

where $a$ represents semi-major axis value and $b$ represents semi-minor axis value.

One seeks to determine the intersection between an arbitrary right circular cone (with parameters corresponding to typical values obtained from a UAV or LEOS) and the WGS84 ellipsoid.

This intersection will in general be a curve, or more accurately, may comprise one or more disjoint curves. To see that for example two (non-intersecting) curves could be obtained, suppose that the apex of the cone lay outside the ellipsoid, say above the north pole, with the axis of the cone passing through the centre of the ellipsoid. Then if the semiangle of the cone were sufficiently small, the cone would intersect the upper hemisphere of the ellipsoid in a closed curve, and separately it would intersect the lower hemisphere of the ellipsoid in a second closed curve. Of course, this is not the only possibility.

For the purposes of establishing and recording the intersection of the cone and the ellipsoid, one could imagine choosing a set of neighboring values of $x$ (or $y$ or $z$), and for each value determining the real intersection points (which will normally be 0, 2 or 4 in number) of the two equations. (We say normally, because the number may be different if there is a point at which the two surfaces touch and have some sort of common tangency.)

Given the possibility of the intersection comprising two separate curves, it would be important to determine, for neighboring values of $x$, which intersection points were on the same curve, and which were not. When the intersection points are remote from one another, this will be easy. It will be less straightforward if curves merge, touch, etc.

IV-A Find the intersections of a cone with an ellipsoid

To parametrize a right circular cone, we can divide the cone surface into two parts, an apex $\mathbf{r}=(r_{x},r_{y},r_{z})$ and a set $\Phi$ of straight rays which start from the apex and ends in the infinity. Suppose $\alpha_{c}$, $\beta_{c}$ and $\gamma_{c}$ are direction cosines of one of the rays in $\Phi$. The ray equation can be given by

[TABLE]

where $s\in[0,\infty)$ is the range of the ray. The constraint of non-negative value on $s$ ensures that (41) describes a right circular cone rather than a double cone.

Now the problem of finding the intersections of two quadratic surfaces (a right circular cone and an ellipsoid) has been transformed to finding the intersections of a set of rays described by (41) with the ellipsoid (40). According to (41), it is observed that the direction cosines $\alpha_{c}$, $\beta_{c}$ and $\gamma_{c}$ can be obtained by transforming the direction cosines of the rays on the surface of the original cone (29). We first pick one ray on the surface of the cone (29) and identify its direction cosines $\alpha_{r}$, $\beta_{r}$ and $\gamma_{r}$ as follows:

[TABLE]

where $\zeta$ denotes the angle from the ray to its projection on $x-y$ plane, and $\eta$ denotes the angle between the $x-{\text{axis}}$ and the projection of the ray on $x-y$ plane. See Fig. 4 for the illustrations of $\zeta$ and $\eta$. Let $\mathbf{d}_{r}:=[\alpha_{r}~{}\beta_{r}~{}\gamma_{r}]^{\top}$ be an unit direction vector of an arbitrary ray on cone (29). The unit direction vector $\mathbf{d}_{r}:=[\alpha_{t}~{}\beta_{t}~{}\gamma_{t}]^{\top}$ of the mapping of the ray on the transformed cone $\eqref{eq:transformed_cone}$ can be directly calculated by

[TABLE]

Here $\mathbf{R}$ is the rotation matrix defined in (31).

Then for any particular ray, $r_{x},r_{y},r_{z},\alpha_{c},\beta_{c}$ and $\gamma_{c}$ are known. By substituting (41) into the ellipsoid equation (40), we have a second-order equation in $s$ with the form

[TABLE]

where

[TABLE]

are all known parameters with $q_{0}=-1$, $q_{1}=q_{2}=1/a^{2}$ and $q_{3}=1/b^{2}$. The real, non-negative solutions of (44) define points on the ray corresponding to its intersections with ellipsoid. Note that the number of solutions of (44) can be zero, one or two, depends on the value of the ”quadratic discriminant”:

[TABLE]

The number of solutions of (44) also represents the number the intersection points of the ray with the ellipsoid.

IV-B Find the intersections of the FDOA cone with the earth’s terrain

The aim of this subsection is to extend the calculation of the previous subsection involving the WGS84 ellipsoid to provide an algorithm to find the intersections of the FDOA cone with the earth terrain and plot these intersection points on a 3D map tool (e.g. google earth). We assume that the vehicle’s measurements of its position, velocity and the Doppler shift are all accurate, i.e. without error.

Since the earth terrain is uneven, we cannot find equations to describe the terrain, which means we cannot directly calculate the intersections of a cone with the terrain by solving mathematical equations. Motivated by the observations that the earth terrain can be modelled by the discrete point set DTED, our idea is to find the intersection of the cone with the WGS84 ellipsoid and map these intersection points to the DTED dataset, according to the rays generated between the cone apex and the intersection points on the ellipsoid surface.

Our algorithm is consisted of the following steps:

Introduce $\mathrm{D}$ as a set of points comprising the DTED. Thus

[TABLE]

Note that $\phi_{i}$ and $\lambda_{i}$ are defined with respect to the WGS84 ellipsoid, while $H_{i}$ is defined with respect to the EGM96 geoid [13] . For further use, we need to transfer $H_{i}$ (orthometric height) to $h_{i}$ (ellipsoid height) by using (25). We denote this new set as $\hat{\mathrm{D}}$. 2. 2.

By applying $\eqref{eq:geotocartesian}$ on each point in $\hat{\mathrm{D}}$, we obtain a new set $\hat{D}_{c}$, in which the coordinates for each point are defined with respect to the ECEF Cartesian coordinate system. 3. 3.

By using the method introduced in the previous subsection, the vehicle determines a set of the intersection points of the FDOA cone (associated with one single measurement) with the WGS84 ellipsoid. Introduce

[TABLE]

to denote the set of intersection points. Note that the coordinates associated with the points in $\mathrm{A}$ are defined with respect to the ECEF Cartesian coordinate system. 4. 4.

Select a point $p_{\mathrm{A}}^{i}$ in $\mathrm{A}$. By associating the coordinates of $p_{\mathrm{A}}^{i}$ with the receiver coordinates $\mathbf{r}$, we can define a line described by the following equation

[TABLE]

where $a_{i},b_{i},c_{i}$ can be obtained by substituting the coordinates of $p_{A}^{i}$ and $\mathrm{r}$ into the above equation. For later use, we denote this line as $l_{i}$. It is one of the rays that comprise the surface of the FDOA cone. 5. 5.

Now we have a transformed DTED point set $\hat{\mathrm{D}}_{c}$ and the line $l_{i}$. For each point $p_{k},k=1,2,\ldots$ in $\hat{\mathrm{D}}_{c}$, we can calculate its Euclidean distance $d_{k}^{i}$ to $l_{i}$. We then select a fixed threshold $tr_{i}$ (the threshold is determined according to the spacing accuracy of the DTED sets) and figure out a new point set $\hat{D}_{i}$, defined as

[TABLE]

The mapping of $p_{A}^{i}$ onto $\hat{D}_{c}$ is thus included in $\hat{D}_{i}$. 6. 6.

For each point $\hat{D}_{i}$, we calculate its Euclidean distance to the receiver, the point $\tilde{p}$ associated with the minimum distance to the receiver is the mapping point of $p_{A}^{i}$ onto $\hat{\mathrm{D}}_{c}$. 7. 7.

By repeating the steps 4)-6) for each $l_{i}$, we can finally find the cone-terrain intersection points. 8. 8.

For each cone-terrain intersection point, we use (3) to transform its ECEF coordinates into the geographical coordinates to depict it in a 3d map tool.

In the following, we use an example associated with UAV cases to illustrate the performance of our algorithms. The UAV position is selected to be close to the city of Adelaide, Australia. The FDOA cone, the intersections of the FDOA cone with the earth terrain, and the UAV are depicted in Fig. 5. We also provide Fig. 6 for the purpose of clear observations on how the intersection points are distributed on the earth’s terrain.

V Error analysis

This section aims to discuss the curve shift associated with different errors raised in practical scenarios. We identify two potential error sources, 1) the fact that the nominal frequency of the emitter may be different from its true frequency; 2) the changing refractive index of the atmosphere, which results in a variation of the speed of the electromagnetic wave. We provide numerical examples to illustrate the curve shifts associated with different errors. Since UAV and LEOS are operated in different environments, we divide the discussions into two cases. For LEOS case, we further consider the relativistic Doppler effect. For all the examples in this section, we assume that the measured frequency at the receiver and the navigation (position, velocity and attitude) data of the vehicle are accurate, i.e. without errors. Note that for calculation simplicity, we assume that the nominal frequency of the emitter is 299,792,458Hz.

V-A Error sources

We now discuss the first error source. The carrier frequency of the emitter is generated according to its built-in oscillators. The output frequency of an oscillator has a drift due to its internal changes with respect to time plus changes in the environment. So there is potentially an offset between the nominal frequency and the true frequency of the emitter. Before this section, we assumed that the true frequency rather than the nominal frequency of the emitter is precisely known to the receiver. This can only be assumed if, for example, there is a separate stationary receiver providing the information to the moving receiver, or the emitter is actually reflecting signals due to its being illuminated by a radar, with the position, velocity and precise frequency of the radar transmitter known to the moving receiver.

We then explain the second error source. It is well-known that the speed of electromagnetic wave in a medium different to vacuum, is different to its speed in the vacuum. Since the speed of electromagnetic wave changes with respect to the air refractive index when travelling in the air, a small error will be introduced when calculating the semi-angle of the FDOA cone using (26). The speed of electromagnetic wave in the air is normally calculated using:

[TABLE]

where $n$ is the refractive index of the air. The refractive index of the air varies according to air temperature and pressure, which leads to a refractive index gradient in the atmosphere [15]. Note that air temperature and pressure can also be described with respect to the altitude. In this paper, we refer to an air refractive index model described by a function of altitude proposed in [16].

In the following, we will use numerical examples to illustrate the curve shift associated with different errors. In all the figures, the curves that represent the intersections of the true FDOA cones with the WGS84 ellipsooid are marked in yellow and the curves calculated with respect to the errors are marked in red.

V-B The UAV cases

According to the U.S. Department of Defense [17], UAVs can be conveniently classified into five categories according to their sizes (see Appendix for the classifications) . In our application scenarios, a medium- or large-size UAV is capable of carrying on the tasks. In this subsection, we assume that the UAV speed and height are 50m/s (180km/h) and 2000m, respectively. The yaw angle of the UAV is selected as [math], the pitch angle is selected as [math] and the roll angle of the UAV is selected as [math]. We consider three situations associated with different magnitudes (close to [math], close to [math] and close to [math]) of the semi-angle of the true FDOA cone, to better describe the situations in practical scenarios.

V-B1 Curve shift associated with the error between the nominal frequency and the true frequency

• •

UAV example 1: Assume that the measured frequency at the receiver is 299,792,501.33Hz, and the true frequency of the emitter is 299,792,468Hz. There is a 10Hz error between the nominal frequency and the true frequency of the emitter. The semi-angle of the true FDOA cone is calculated according to the true frequency of the emitter, which can be obtained as [math]. Since the receiver has only the knowledge of the nominal frequency of the emitter, the semi-angle of the FDOA cone calculated at the receiver is [math]. We plot the two intersection curves in Google Earth, see Fig. 7, and observe that the magnitude of the curve shift ranges from approximately 730m to 4,000km 333Note that we do not consider the fact that the signal from the emitter can be blocked due to the eccentricity of the earth surface. The curve shift does not represent the error of the position of the emitter..

• •

UAV example 2:

We still assume that the true frequency of the emitter is 299,792,468Hz. By assuming that the measured frequency at the receiver is 299,792,507.93Hz, the Doppler shift calculated according to the nominal frequency and the true frequency of the emitter are 49.93Hz and 39.93Hz, respectively. The semi-angle of the FDOA cone associated with the true Doppler shift is [math]. The semi-angle calculated according to the difference between the measured frequency and the nominal frequency is [math]. The two intersection curves for this example are depicted in Fig. 8. The minimum magnitude of the curve shift is approximately 2,250m and the maximum magnitude of the curve shift is too large to be determined in Google Earth. This is because the size of the yellow curve is very small (the radius is about 800m) and the red curve in Fig. 8 has a similar shape as the red curve in Fig. 7, with a much bigger size.

• •

UAV example 3:

Assume that the measured frequency at the receiver is 299,792,470.615Hz, and the true frequency of the emitter is 299,792,468Hz. The true Doppler shift is 2.615Hz, thus the semi-angle of the true FDOA cone is [math]. The Doppler shift calculated according to the nominal frequency of the emitter is 12.615 Hz and its associated FDOA cone has a semi-angle with the value of [math]. The magnitude of the curve shift measured in Goolge Earth ranges from approximately 500m to 2,600km

By observing the above three UAV examples, we can conclude that there is a significant curve shift associated with the error between the nominal frequency and the true frequency of the emitter. When applying our algorithm on the UAV platforms, one must be very careful to ensure that the knowledge of the transmission frequency of the emitter is accurate enough.

V-B2 Curve shift associated with the changing refractive index of the atmosphere

The UAV operation altitude in our application scenarios is less than 5km. According to the air refractive model proposed in [16], the air refractive index can be regarded as a constant with a value 1.0003. We will use this value in the following examples.

• •

UAV example 4:

Assume that the Doppler shift measured at the receiver is 43.3Hz. Then the semi-angle of the FDOA cone calculated using the speed of electromagnetic wave in the vacuum is [math]. By considering the refractive index of the air, we can directly substitute (49) into equation (26) and obtain that the semi-angle of the FDOA cone considering the refractive index of the air is [math]. We depeict the two intersections curves in Fig. 10. The range of the curve shift varies from 3m to 6,500m.

• •

UAV example 5:

In this example, we assume that the Doppler shift measured at the receiver is 49.93Hz. Thus the semi-angle of the FDOA cone calculated using the speed of the electromagnetic wave in the vacuum is [math]. Considering the speed of the electromagnetic wave in the air, the new magnitude of the semi-angle of the FDOA cone can be obtained as [math]. Fig. 11 shows the two intersection curves. The curve shift of the two curves varies from 20m to 40m.

• •

UAV example 6:

In this example, we assume that the Doppler shift measured at the receiver is 2.615Hz. The semi-angle of the FDOA cone calculated using the speed of an electromagnetic wave in a vacuum is [math]. Considering the speed of an electromagnetic wave in the air, the magnitude of the semi-angle of the true FDOA cone is [math]. The curve shift of the two curves varies from 3m to 1,500m.

According to the three examples for UAV cases considering the changing refractive index of the atmosphere, we can find that the magnitude of the curve shift is actually minor. Though the maximum curve shifts in example 4 and 6 are 1500m and 3500m, they both appear at the positions far away from the receiver (the distance can be more than 4,000km). And it is obvious that signal emitted at those positions cannot be received by the receiver. In conclusion, we may neglect the effect of the changing refractive index of the atmosphere when applying our algorithm in practice.

V-C The LEOS cases

V-C1 Curve shift associated with the error between the nominal frequency and the true frequency of the emitter

In all the examples for LEOS, we assume that the LEOS is above the equator and the roll, pitch and yaw are all selected as [math]. We also assume that the satellite speed is 7800m/s (28,080km/h). The height of the satellite is 200,000m (200km). These two parameters are typical for LEOS. In all the examples provided below, we assume that the error between the nominal frequency and the true frequency is 60Hz.

• •

LEOS example 1: The true transmission frequency of the emitter is 299,792,518Hz. Suppose that the measured frequency at LEOS is 299,798,033.4329Hz. Thus the true Doppler shift is 5515.4329Hz. The semi-angle of the FDOA cone calculated according to the true frequency of the emitter is [math]. The semi-angle of the FDOA cone calculated according to the nominal frequency of the emitter $45.62^{o}$. The two curves are shown in Fig. 13. The curve shift of the two curves varies from approximately 5,000m to 15,000m.

• •

LEOS example 2:

Since the LEOS orbit is very close to circular and the moving direction of LEOS is close to be parallel to the plane of the horizon. If we consider the situation that the semi-angle of the FDOA cone is close to [math], then it is impossible to find an emitter on the earth’s surface to generate such the FDOA cone. We do not need to provide illustrations to this example.

• •

LEOS example 3: We still assume that the true transmission frequency of the emitter is 299,792,518Hz. Assume that the Doppler shift measured at the LEOS is 347.94Hz. Because of the error between the nominal frequency and the true frequency, the output of the Doppler shift at the LEOS is 407.94Hz. Then the semi-angle of the FDOA cone calculated according to the nominal frequency of the emitter is [math]. The semi-angle of the FDOA cone calculated according to the true transmission frequency of the emitter is [math]. The two curves are depicted in Fig. 12. The curve shift of the two curves varies from 2,000m to 100,000m.

By observing the above LEOS examples, we can find out that the magnitude of the curve shift is not negligible, though the change of the semi-angle of the FDOA cone is minor. This is due to the fact that LEOS are operated at a high altitude. When applying our algorithms on the LEOS, it is still required to have a accurate knowledge of the true transmission frequency of the emitter.

V-C2 Curve shift associated with the changing refractive index of the atmosphere

The refractive index of air varies according to a function of temperature and pressure, which leads to a refractive index gradient in the atmosphere. Although this gradient is very small, atmospheric refraction of electromagnetic waves is observable due to the large distances traveled in the atmosphere between the transmitter and the emitter.

The earth atmosphere can be divided into several layers, including Troposphere, Stratosphere, Mesosphere and Ionosphere. From the viewpoint of atmospheric refraction, only the first two layers, the troposphere and the stratosphere, are important. In the upper layers of the atmosphere the air is so rarefied that the refractive index can be considered to be unity.

In this subsection, we roughly estimate the curve shift by assuming a simple two layer atmosphere model instead of the complicated model described in [16]. Our model is described as follows: assume that the height of the Troposphere is $20km$, and the Stratosphere starts at $20km$ and ends at $50km$. The refractive index is assumed to be a constant for each layer, i.e, the refractive index $n_{s}$ of the Stratosphere is equal to 1 and the refractive index $n_{t}$ of Troposphere is equal to 1.0003.

Now we describe our scenario to calculate an estimated curve shift. The graphical illustration of this scenario is in Fig. 15. Suppose the LEOS is moving horizontally and its calculated semi-angle of the FDOA cone associated with one single measurement is [math]. According to the above 2-layer atmosphere model and the well-known Snell’s law, it is straightforward to obtain that the curve shift is about 41m, which is small when compared to the curve shift associated with the error between the nominal and the true frequencies.

V-C3 Curve shift associated with the relativistic Doppler effect

The relativistic Doppler effect is the change in frequency of electromagnetic waves, caused by the relative motion of the source and the observer (as in the classical Doppler effect), when taking into account effects described by the special theory of relativity. The relativistic Doppler effect is different from the non-relativistic Doppler effect as the equations include the time dilation effect of special relativity and do not involve the medium of propagation as a reference point.

Due to relativistic effects, clocks on the receiver are time dilated relative to clocks at the source. The measured frequency in the receiver, denoted by $f_{r}$, can be described by

[TABLE]

where $f$ is the classical Doppler shift described by (27). $\rho$ is called the Lorentz factor, which is defined as

[TABLE]

in which $\varsigma:=v_{r}/c$ is the velocity of the receiver in terms of the speed of light.

We now use an example to illustrate the curve shift associated with the relativistic Doppler shift. We use the parameters selected in LEOS example 1 in this example. We assume that there exists a FODA cone with an accurate semi-angle of [math] associated with a single Doppler shift measurement at the LEO satellite. Now we consider the relativistic Doppler effect. According to (50), it can be obtained that $\rho$ is equal to a value that is very close to 1 with the error less than $10^{-10}$. By using this value, we can obtain that the change of the semi-angle of FDOA cone is $1.94\times 10^{-7}$ degrees. By assuming that the earth surface is flat and there is no refraction through the air, it can be directly calculated that the curve shift is about $0.0001352m$, which is small enough to be neglected.

VI Conclusions and future work

This paper investigates the problem of determining a curve on the earth’s surface on which a stationary emitter must lie based only on the measured frequency on a mobile vehicle (UAV or LEOS) and the vehicle’s navigation data (position and velocity). Our investigated scenario considers the following two facts: 1) Coriolis effect; 2) the bumpy earth’s surface. We prove that the relative velocity component due to the Coriolis effect makes zero contribution to the Doppler shift. We then provide methods for building equations to describe a FDOA cone and finding its intersection curve with the WGS84 ellipsoid. The curve comprises points, which are mapped into a DTED dataset to find the curve on the earth’s terrain. We further provide numerical examples to illustrate how the errors resulting from the nonconstant refractive index of the atmosphere and from lack of precise knowledge of the transmitter frequency affect the positions of curves. Future work will focus on estimating the true frequency of the emitter using multiple frequency measurements and their associated intersection curves on the surface of the WGS84 ellipsoid.

-A Abbreviations

[TABLE]

-B Notations

[TABLE]

-C UAVs Classification

[TABLE]