Machine and process for self localization using doppler

Abstract

A system, method and apparatus is provided for calculating the position of a single receiver by receiving a signal from an emitter. Either the receiver or the emitter may be stationary relative to the surface of the Earth. The invention may use a particle filter to optimize the efficiency of the position calculation. The invention may further utilize a particle filter to estimate the frequency drift and bias of the emitter and/or receiver, and use the estimated frequency drift to improve the position calculation. In one embodiment, the invention may be used to calculate the position of a single GPS receiver given the position of, and a transmission from, a single GPS satellite.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

The present application claims benefit of similarly titled U.S. Provisional Patent application Ser. No. 61/489,953, filed on May 25, 2011, which is incorporated by reference in its entirety. The present application also claims the benefit of U.S. Provisional Patent application Ser. No. 61/519,597, titled “Diversified Doppler for Single Platform Geolocation,” filed on May 25, 2011, which is incorporated by reference in its entirety. This application is related to contemporaneously filed U.S. patent application Ser. No. 13/479,573, titled “Diversified Doppler for Single Platform Geolocation,” filed May 24, 2012, which is incorporated by reference in its entirety.

CROSS-REFERENCE TO RELATED TECHNICAL DISCLOSURES

The present application makes use of technology disclosed in these related references, each of which is hereby incorporated by reference in its entirety:

• • Witzgall, Hanna—A Reliable Doppler Based Solution for Single Sensor Geolocation, Submitted to IEEE Military Communications Conference (MILCOM) for publication in October 2012.
  • Witzgall, Hanna, et. al.—Single Platform Passive Doppler Geolocation with Unknown Emitter Frequency, Aerospace Conference, March 2010 IEEE, ISBN 978-1-4244-3887-7.
  • Witzgall, Hanna, et. al.—Doppler Geolocation with Drifting Carrier, Military Communications Conference, November 2011 IEEE, ISBN 978-1-4673-0079-7.
  • Arulampalam, M. Sanjeev, et al.—A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking, IEEE Transactions on Signal Processing, February 2002, ISSN 1053-587X.

FIELD OF THE INVENTION

The present invention relates to the localization of a radio frequency (“RF”) receiver. More particularly, the present invention relates to self localization of an RF receiver using Doppler information collected from one or more moving beacons. If the receiving platform's local oscillator contains a changing measurement error the localization solution becomes difficult to solve. This invention describes a technique that addresses this problem and allows inexpensive receivers to accurately self localize from one or more moving beacons.

BACKGROUND

Frequency of Arrival (FOA) measurements provide bearing angles that can be used for localization provided the known position, velocity and frequency of the moving emitter platform. In the well-known Doppler equation, the frequency observed by a receiving platform is shifted in a nonlinear fashion relative to the emitter's carrier frequency f_e and the geometric relationship between the receiver and the emitting RF platform, is shown below in Equation 1:

$\begin{array}{cc}\Delta f=\left(\frac{v\cos \theta }{c}\right)f_e& \mathrm{Equation}1\end{array}$
where Δf the Doppler shift, v is the speed of the receiver relative to the emitter platform, θ is the cone angle between the receiver and the platform velocity vector, and c is the speed of light. Multiple Doppler observations can be triangulated and used to solve for the unknown receiver location.

To date, the most common methods for geolocation employ multilateration (i.e., multiple emitters). Typically these methods employ measurements of Time Difference of Arrival (“TDOA”) and/or Frequency Difference of Arrival (“FDOA”) of RF signals received from at least two emitters. For example, the 3D position of a receiver may be determined by measuring the TDOA of a signal from four or more spatially separate emitters located different distances from the receiver. If the emitting platform and receiver are in relative motion with respect to one another, FDOA may be employed either independently or in addition to TDOA to determine the receiver position based on \ the observed Doppler information. The popular global positioning system (GPS) uses the TDOA measurements collected from four or more satellites to act as RF beacons for a GPS receiver.

FDOA uses the difference of frequencies collected simultaneously from at least two different moving emitters. As shown in Equation 2, below, the FDOA observable Δω describes the Doppler shift caused by the relative change in range between the platform and two emitters:

$\begin{array}{cc}\Delta \omega =\frac{f_e}{c}\frac{d}{dt}\left(r_1-r_2\right)& \mathrm{Equation}2\end{array}$
where the variable r_k represents the range between the receiver k and an emitter.

A drawback of FDOA-based receiver location algorithms is that they require at least two or more simultaneous emitters. Furthermore, local oscillators (sensors) used in receivers such as cellular phones or handheld GPS devices often impart a bias and drift into their frequency measurements. For example, a given sensor may have an RF measurement bias, meaning an offset between the reported measurement and the actual frequency sampled. Further, this offset may have a stationary component, and also a component that drifts or changes over time. For example, as the temperature of the sensor changes, the offset may increase or decrease. This can cause a large localization error if not corrected.

This correction is sometimes done using known calibration tones present in the collected data, however, this requires even more resources, including a potentially cooperative emitter, and further motivates a robust, single beacon, single receiver FOA localization solution.

The FOA measurement equation is derived directly from Equation 1 by rearranging terms to isolate the observed Doppler-shifted frequency f_obs from the unobserved Doppler shift Δf that depends on the emitter's carrier frequency f_e. Furthermore, the cosine term can be expanded to explicitly contain the state variables of receiver location [x_r,y_r]^T and emitter frequency f_e by writing out the dot product between the emitter velocity [v_x,v_y]^T and the line of sight between the receiver and the emitting platform [x,y]^T:

$\begin{array}{cc}{f}_{\mathrm{obs}}=\left[1-\frac{1}{c}\frac{\left(x-x_r\right)v_x+\left(y-y_r\right)v_y}{c\sqrt{{\left(x-x_r\right)}^2+{\left(y-y_r\right)}^2}}\right]f_e-f_r& \mathrm{Equation}3\end{array}$

A primary reason that the FOA method is generally not considered an operational localization option is due to its tendency to produce unreliable results.

One reason for this is the fact that many commercial receiver systems do not have strict tolerances on oscillator stability which means that a receiver's frequency estimate can be both biased and drift significantly during normal operation. Failing to account for frequency drift observed over the course of triangulation time required to determine a unique solution will cause gross inaccuracies in location estimates. However, including a non-stationary drift term in the location solution eliminates the suitability of Gauss Newton and similar techniques, because it requires a stationary state model. This motivates the need for a new solution that can efficiently represent and track in a non-linear, non-convex, multimodal state space.

SUMMARY OF THE INVENTION

Exemplary embodiments described herein overcome the above described drawbacks of conventional multi-emitter FDOA localization techniques. In certain embodiments, the described invention allows for rapid localization of a single receiver using a single moving emitter platform. These embodiments overcome the location inaccuracies of conventional FOA methods by solving simultaneously for the receiver location and any potential frequency drift imparted by either the receiver in its measurements or the known beacon's transmission frequency. Specific embodiments may utilize particle filtering algorithms to recursively update multimodal state densities that are typical of solutions involving both unknown receiver location and non-stationary receiver measurement frequency drift. These embodiments provide a more robust and accurate solution compared with a conventional implementation of a FOA method.

In certain embodiments, multiple independent observables based on Doppler characteristics may be utilized to provide instantaneous or near-instantaneous estimates of a receiver location, velocity, and/or carrier frequency measurement bias. These measurements may be employed alone or in combination with particle filtering algorithms to achieve unprecedented accuracy and speed using single-emitter, single-platform geolocation.

In one aspect of the invention, a method of determining a location of a radio frequency (RF) receiver using an RF receiver located on a single platform, and an emitter in motion relative to the receiver, is provided. The method includes receiving, by the receiver, an RF signal transmitted by the emitter over a period of time. Frequency of arrival (FOA) data of the RF signal, such as but not limited to Doppler-shifted frequency of the emitter, a location of the emitter, and a velocity of the emitter, may be measured, received or known over the period of time. This FOA data may be employed to determine an estimate of a frequency measurement bias of the receiver, wherein the estimated frequency measurement bias of the receiver includes an estimate of a frequency drift. Finally, an updated estimate of the location of the receiver may be determined based on the estimate of the frequency measurement bias of the receiver.

In another aspect of the invention, a system capable of determining the location of a receiver is provided. The system includes a platform containing a receiver, and a processor in communication with the receiver. The receiver is adapted to receive an RF signal transmitted by an emitter over a period of time. The processor is adapted to measure frequency of arrival (FOA) data (e.g. Doppler-shifted frequency of the emitter, a location of the emitter, and a velocity of the emitter) of the RF signal over the period of time. The processor may use this data to determine an estimate of a frequency measurement bias of the receiver based on an estimate of the location of the receiver, wherein the estimated frequency measurement bias of the receiver includes an estimate of a frequency drift. Finally, the processor may determine an updated estimate of the location of the receiver based on the estimate of the frequency measurement bias of the receiver.

In yet another aspect of the invention, a platform capable of self localization is provided. The platform may include a receiver for receiving an RF signal from an RF emitter. The platform may also include a processor for determining an estimated position of the RF receiver based on two or more measurements of the received RF signal selected from the group consisting of FOA, FOA-DOT, FOA-DOT-DOT, FDOA, and D-FDOA.

These and other aspects of the invention will be better understood by reading the following detailed description and appended claims.

BRIEF DESCRIPTION OF THE FIGURES

The preferred embodiments of the present invention are illustrated by way of example and are not limited to the following figures:

FIG. 1 illustrates a system according to one exemplary embodiment of the invention;

FIG. 2 illustrates an exemplary single emitter geolocation method according to an embodiment of the invention;

FIG. 3 depicts a graphical representation of multiple carrier frequency estimates over a platform flight path time;

FIG. 4 illustrates a system according to one exemplary embodiment of the invention;

FIG. 5 illustrates an exemplary flight path of a platform with respect to the location of an emitter;

FIG. 6 illustrates the carrier frequency drift of the emitter of FIG. 5.

FIG. 7 depicts a graphical comparison of an observed Doppler effect for a stable emitter 701 to the observed Doppler effect 702 of an emitter having the carrier frequency drift of FIG. 6;

FIG. 8 depicts the geolocation accuracy of 10 test runs of a system running an exemplary DPF algorithm compared to the FOA Cramer-Rao lower bound (CRLB), assuming an unknown but stationary carrier;

FIG. 9 depicts the error of an exemplary DPF algorithm and a Gauss Newton algorithm on a log 10 scale;

FIG. 10 illustrates an exemplary flight path of a platform with respect to the location of an emitter;

FIG. 11 depicts a graphical comparison of an observed Doppler effect for a stable emitter 1101 to the observed Doppler effect 1102 of an emitter having the carrier frequency drift of FIG. 6;

FIG. 12 depicts the error of an exemplary DPF algorithm and a Gauss Newton algorithm on a log 10 scale; and

FIG. 13 depicts a graphical comparison of FDOA CRLB using multiple platforms to a FOA CRLB employing a single platform.

DETAILED DESCRIPTION

Various embodiments and aspects of the invention will be described with reference to details discussed below, and the accompanying drawings will illustrate the various embodiments. The following description and drawings are illustrative of the invention and are not to be construed as limiting the invention. Numerous specific details are described to provide a thorough understanding of various embodiments of the present invention. However, in certain instances, well-known or conventional details are not described in order to provide a concise discussion of embodiments of the present inventions. All terms used herein are intended to have their ordinary meaning in the art unless otherwise provided.

The following table contains a list of acronyms and abbreviations with corresponding meaning that are discussed throughout this disclosure.

TABLE 1

Acronym     Meaning

AOA         Angle of Arrival
CAF         Cross-Ambiguity Function
CRLB        Cramer-Rao Lower Bound
CONOP       Constrained Optimization
D-DFOA      Difference-DFOA
DFOA        Difference of FOA
FDOA        Frequency Difference of Arrival
FOA         Frequency of Arrival
FOA-DOT     Rate of Change of FOA
FOA-DOT-DOT Doppler Curvature = 2nd Derivative of
            FOA
TDOA        Time Difference of Arrival

Referring to FIG. 1, an exemplary system according to the invention is illustrated. As shown, a receiving platform 100 is on the surface of a planet 101. An emitter 110/112 is moving with respect to the platform 100. The platform 100 may be any equipment capable of including a suitable receiver and accompanying hardware, firmware and software, while moving with respect to an emitter 110/112, but may be either moving or stationary with respect to the surface of the planet 101. For example, a platform 100 may include, without limitation, a car, truck, van, bus, military vehicle, plane, UAV, helicopter, boat, bike, motorcycle, scooter, train, cell phone, GPS receiver, cruise missile, self-guided munition or the like. It will be appreciated that a platform 100 may also refer to a wearable system, such as a backpack or the like worn by a person. In preferred embodiments, the platform 100 may be an aerial vehicle such as a UAV or plane, or a handheld device such as a cell phone.

The platform 100 may comprise a receiver for receiving RF signals from the emitter 110/112. In certain embodiments, the receiver may comprise one or more antennas (e.g. omnidirectional antenna) adapted to receive RF signals from an emitter 110/112 of interest and may be employed to determine one or more Doppler measurements for the emitter (e.g., FOA). In certain embodiments, the platform 100 may also include a GPS receiver for providing time and position reference information.

The platform 100 will typically include an on-board computer in communication with the receiver. The computer comprises at least a processor for processing data and a memory for storing data, and may be adapted to process the received emitter signals in order to generate an estimate of the platform's 100 location. The received signals, estimated emitter location and/or a platform location may be stored in volatile or non-volatile memory, or in the registers of the computer, and such information may be updated in real time or near real time.

In certain embodiments, the computer may be in communication with a ground-based server (not shown) via a network or the like such that data may be transmitted from the computer of the platform to the server. The server and/or the computer may comprise a processor for processing data, memory, I/O ports for receiving data, and may be coupled to any number of user input or output devices such as a keyboard, mouse, monitor, printer or the like. The received signals, emitter or platform location and/or estimated emitter or platform location may be displayed on, for example, a monitor or other display in communication with the computer and/or the server. Although only a single platform 100 is shown, a server may be in communication with any number of computers detecting any number of emitters 110/112.

In certain embodiments, the computer or server may be in communication with an emitter database (not shown), which may be located either at the computer as, for example, a local or external drive, or at the server or other remote location in communication with the computer. Generally, the emitter database may refer to a database used to store information about emitters 110/112 tracked by a platform's computer. If multiple platforms 100 are tracking multiple emitters 110/112, the emitter database may comprise emitter information collected from each of the platforms.

In operation, the platform 100 receives RF signals from one or more emitters 110/112 located at a known or unknown location. Examples of emitters include, but are not limited to cellular phones (e.g., smart phones), GPS devices, GPS or other telemetry satellites, radios, audio/video equipment, personal communications systems, and military and commercial radar and communication equipment.

In certain embodiments, the above described system may employ a particle filtering algorithm to recursively update multimodal state densities related to platform or emitter location and/or emitter carrier frequency. For example, the inventive methods may solve the complex multimodal state equations that arise with a drifting RF carrier for the RF receiver location. Particle filters are sequential Monte Carlo methods based on point mass (or “particle”) representations of probability densities, which can be applied to a state-space model and which may generalize a Kalman filtering method. However, the invention is not necessarily limited to using particle filters to solve for the multimodal state densities found in the simultaneous solution of receiver location and frequency measurement bias estimate, and non-linear optimization techniques can also be employed to fit the data in a piece-wise fashion. For completeness, a particle implementation will be described.

Particle filters are generally concerned with representing the posterior probability density p(x_1:k|z_1:k) using a set of discrete particles x_kⁱ and their associated weights w_kⁱ, as shown below in Equation 4:

$\begin{array}{cc}p\left(x_{1:k}|z_{1:k}\right)\approx \sum _{i=1}^N{w}_k^i\delta \left(x_k-x_k^i\right)& \mathrm{Equation}4\end{array}$
where δ is the Dirac function, the particle's weight w_kⁱ represents its conditional likelihood, and N is the number of particles.

Particles x_kⁱ may be drawn from an importance or proposal density q according to the principle of importance sampling.

$\begin{array}{cc}{x}_k^i\sim q\left(x_k^i|x_{1:k-1}^i,z_k\right)& \mathrm{Equation}5\end{array}$

As shown in Equation 5, importance sampling acknowledges the difficulty of directly sampling from the often complex posterior probability density p(x_0:k|z_1:k) and instead proposes to approximate the true posterior with an importance or proposal density q(x_0:k|z_1:k) from which it is relatively simple to sample (e.g., assuming a Gaussian distribution).

Weights may be employed to correct for the difference in probability between the samples drawn from the proposal density and the true posterior density. For example, a particle filter may normalize updated weights and may send these weights to a resampling algorithm. Such resampling of particles helps to avoid the degeneracy condition where only a single particle has any likelihood probability. Particles may be selected with replacement from the probability distribution of the weights' cumulative sum. By eliminating low-probability particles and duplicating high-probability particles, a particle filter can re-distribute particles to relevant regions and improve accuracy with little increase in computation.

Accordingly, a particle filter may comprise four parts. The first part may calculate a proposal density; the second part may sample this proposal density; the third step may calculate the correction weights; and the fourth step may resample. The procedure may be repeated with the resampled data. In this way, a particle filter can provide a recursive approximation of the state density conditioned on the sequence of observations. The maximum a posteriori estimate may be obtained from the approximated posteriori probability by either computing the mean or by taking the sample with the highest weight.

The inventive systems and methods described herein exploit certain properties of particle filters to provide a geolocation solution given a complex multimodal state space composed of receiver location and a non-stationary frequency measurement bias required for FOA. In certain embodiments, solutions are proposed based on the factorization of the conditional joint state probability of the receiver p(xr_kⁱ, fr_kⁱ|z_k) into a location term p(xr_kⁱ|z_k) and a frequency term p(fr_kⁱ|xr_kⁱ,z_k) using Bayes rule as shown in Equation 6:

$\begin{array}{cc}p\left({\mathrm{xr}}_k^i,{\mathrm{fr}}_k^i|z_k\right)=p\left({\mathrm{fr}}_k^i,|{\mathrm{xr}}_k^i,z_k\right)p\left({\mathrm{xr}}_k^i|z_k\right)& \mathrm{Equation}6\end{array}$
where xr_kⁱ=[x_rⁱ, y_rⁱ]^T represents the receiver location vector and fr_kⁱ=[fr_kⁱ, {dot over (f)}r_kⁱ]^T is composed of the current frequency estimate and its bias and drift. Importantly, the joint density cannot be simply factored as shown in Equation 7, because of the conditional dependences of xr_kⁱ and fr_kⁱ.

$\begin{array}{cc}p\left({\mathrm{xr}}_k^i,{\mathrm{fr}}_k^i|z_k\right)=p\left({\mathrm{xr}}_k^i|z_k\right)p\left({\mathrm{fr}}_k^i|z_k\right)& \mathrm{Equation}7\end{array}$

One advantage of the factorization of Equation 6 is that it allows the particles to be drawn from a proposal density approximating only the receiver location density, rather than the joint space also containing receiver frequency measurement bias. The particle's frequency information may therefore be conditionally updated (e.g., using a Kalman filter) based on the receiver location and observation. The representation of the state probability density into a factorized portion based on the drawn locations xr_kⁱ and conditionally updated frequency information fr_kⁱ enables a more efficient state density selection and leads to a better performance with fewer particles. A detailed explanation of particle filters may be found in S. Arulampalam, S. Maskell, N. Gordon, T. Clapp “A Tutorial on Particle Filters for Online Nonlinear/NonGaussian Bayesian Tracking,” IEEE Transactions on Signal Processing, 2002, incorporated herein in its entirety.

Referring to FIG. 2, an exemplary inventive self-localization method employing a Doppler Particle Filter (DPF) algorithm is illustrated. The method begins with state initialization 201 outside of a particle filter function by modeling a region of interest with particles, for example, in a grid like fashion. The grid may be static or dynamic, and the location uncertainty of each particle (i.e., location covariance) may be determined by the initial grid spacing. However, there are many alternative embodiments that may be employed to initialize particle locations, for example, a priori known preferential locations may be utilized. Additionally, particle states may be initialized with frequencies, and particle location may conditioned upon the assigned particle frequency.

As shown in Equation 8, below, the DPF may be initialized by supplying available FOA observations Z_k and the emitter position x₁ and velocity v₁ data covering a time interval. For example, upon receiving FOA observations from a known emitter, an estimated emitter frequency and location may be set, and the receiver's frequency measurement bias may initially be set to zero.

$\begin{array}{cc}\left[{\left\{x_k^i,{\mathrm{Px}}_k^i,\right\}}_{i=1}^N\right]=\mathrm{DPF}\left[{\left\{x_{k-1}^i,{\mathrm{px}}_{k-1}^i,\right\}}_{i=1}^N,Z_k,x_1,v_1\right]& \mathrm{Equation}8\end{array}$

At steps 202, as FOA observations Z_k are collected, the DPF algorithm may employ a nonlinear Bayesian filter, such as but not limited to a Kalman filter, for each particle to determine the most likely frequency track conditioned on that particle's location xⁱ and the FOA observations. The FOA observations may be entered on a sample by sample basis, block basis, segment basis or over an entire data block.

$\begin{array}{cc}\left[f_k^i,{\mathrm{resid}}_k^i\right]=\mathrm{Kalman}\left[f_{1:k-1}^i,P{f}_{k-1}^i,Z_k,x_{k-1}^i,x_1,v_1\right]& \mathrm{Equation}9\end{array}$

A Kalman filter generally assumes that the posterior density at every time step is Gaussian and, hence, parameterized by mean and covariance. A Kalman filter may first update the a priori state information {circumflex over (x)}_k and covariance {circumflex over (P)}_k according to a process model F and the process noise Q, as shown in Equation 10, below:

$\begin{array}{cc}{\hat{x}}_k=F{x}_{k-1}{\hat{P}}_k={\mathrm{FP}}_{k-1}{F}^{\prime }+Q& \mathrm{Equation}10\end{array}$

The Kalman filter may then provide the best posteriori estimate x_k by weighting the a priori prediction error with the Kalman gain K as shown in Equation 11, where R is the measurement noise.

$\begin{array}{cc}{x}_k={\hat{x}}_k+K\left(z_k-{\hat{z}}_k\left({\hat{x}}_k\right)\right)P_k={\hat{P}}_k-K\left(H{\hat{P}}_k{H}^{\prime }=R\right)K^{\prime }K={\hat{P}}_k{H}^{\prime }/H{\hat{P}}_k{H}^{\prime }+R)& \mathrm{Equation}11\end{array}$

In certain embodiments, the Kalman filter may be used to discover the carrier frequency provided the hypothesized receiver location. Therefore, the process equation may be a two-dimensional linear model as shown in Equation 12

$\begin{array}{cc}F=\left[\begin{array}{cc}1& \Delta T\\ 0& 1\end{array}\right],& \mathrm{Equation}12\end{array}$
which may be used to update the prior frequency state x_k-1=[f_r, {dot over (f)}_r]^T for the specified time interval ΔT. The posteriori estimate may use the Kalman gain to weight the prediction error between the FOA observations z_k and the predicted observation {circumflex over (z)}_k provided the receiver location. It will be appreciated that, in other embodiments, a Gauss-Newton method may be employed in addition to or in place of a Kalman filter.

The next step 203 of the particle filter is to determine the importance weights providing the proposal's correction term:

$\begin{array}{cc}{w}_k^i=\frac{p\left(z_k|x_k^i,F_k^i,Z_{k-1}\right)p\left(x_k^i|x_{k-1}^i\right)}{q\left(x_k^i|X_k^i,Z_k\right)}& \mathrm{Equation}13\end{array}$

An importance weight is defined as the ratio between the true updated state probability p(z_k|x_kⁱ)p(x_kⁱ|x_k-1ⁱ) and the probability of drawing that state from the importance density q(x_kⁱ|x_k-1ⁱ, z_k):

$\begin{array}{cc}{w}_k^i=w_{k-1}^i\frac{p\left(z_k|x_k^i\right)p\left(x_k^i|x_{k-1}^i\right)}{q\left(x_k^i|x_{k-1}^i,z_k\right)}& \mathrm{Equation}14\end{array}$

Here the term p(z_k|x_kⁱ) represents the likelihood of the state x_kⁱ, p(x_kⁱ|x_k-1ⁱ) is the probability of the particle's updated state given the process model, and q(x_kⁱ|x_1:k-1ⁱ, z_k) is the probability of having x_kⁱ been drawn from the importance density.

The likelihood term p(z_k|x_kⁱ, F_kⁱ, Z_k-1) may dictate the probability of observing z_k given the receiver location x_kⁱ, its associated carrier track F_kⁱ, and all the available observations Z_k-1. The weight's denominator is simply the proposal density likelihood q(x_kⁱ|X_k-1ⁱ, Z_k) or in other words the Gaussian probability of drawing x_kⁱ given its parent particle's location and covariance. Therefore the weights may be largely determined by how well the frequency track's process model accommodates the FOA observations.

Referring to FIG. 3, Kalman calculated carrier frequency tracks conditioned on particle location and FOA observations (e.g., observed emitter frequency, platform position, and platform velocity over a time period) are shown. Note that FIG. 3 illustrates simulated experimental data tracking an emitter of unknown frequency and drift, with a platform having a known position and velocity. However, given that this experimental data is closely related mathematically to the problem of self-localization, the illustration provided by FIG. 3 is helpful in understanding self-localization. The various carrier paths produced by hypothesizing particles at different locations are illustrated, where line 301 denotes true emitter drift and line 302 denotes the most likely frequency track corresponding to the process model.

As shown, the high curvature of some paths represents the contortions that the carrier drift would need to perform to produce the observed Doppler measurements. It is unlikely that a device would have a sufficiently unstable emitter to produce these types of drifts. A much more likely model is that the carrier frequency drifts slowly over time and changes its drift rate even more slowly. These carrier frequency characteristics can be built into the Kalman filter process model. The more likely the observations support this type of process model, then the more likely that the emitter is at that location. Thus the particle filter may prefer and preferentially weight particles that have a carrier drift that corresponds well with the slow varying process model. It will be appreciated that, in other embodiments, particle weights may be determined as a residual fitting error of a parametric curve fit to observed Doppler measurements found with an optimization technique on a piecewise portion of the measurements.

Referring back to FIG. 2, the particle weights may then be normalized 204 according to Equation 15.

$\begin{array}{cc}{w}_k^i={\tilde{w}}_k^i/\Sigma {\tilde{w}}_k^i& \mathrm{Equation}15\end{array}$

At step 205, the weighted particles may be entered into any resampling function (e.g., according to Equation 16) where higher weights are more likely to be selected multiple times.

$\begin{array}{cc}\left[{\left\{x_k^j,i^j\right\}}_{j=1}^N\right]=\mathrm{Resample}\left[{\left\{x_k^j,w_k^i\right\}}_{j=1}^N\right]& \mathrm{Equation}16\end{array}$
Such resampling prunes unlikely particles, but may cause an inherent degeneracy where all the particles are reassigned as one or two of the most probable samples.

In order to avoid this degeneracy, a child particle may be drawn from a calculated proposal density of the surviving parent particle. Accordingly, the proposal density may be calculated 206 according to Equation 17 for each particle.

$\begin{array}{cc}q\left(x_k^i|X_{k-1}^i,Z_k\right)=N\left(x_k^i|x_{\mathrm{opt}}^i,P{x}_k^i\right),& \mathrm{Equation}17\end{array}$
where (x_optⁱ, Px_kⁱ)=EKF[f_1:kⁱ, z_k, rex_kⁱ, x₁, v₁].

As shown in Equation 17, the proposal density q(x_kⁱ|X_k-1ⁱ, Z_k) is defined as the Gaussian N(x_kⁱ|x_optⁱ, Px_kⁱ) calculated from an extended Kalman filter (EKF) conditioned on the observations and the particle's associated frequency track.

At step 207 a sample child particle x_kⁱ˜N(x_kⁱ; x_optⁱ, Px_kⁱ) may be drawn from the calculated proposal density. In this fashion, the particles are efficiently herded into high probability regions of the state space and degeneracy is largely avoided.

It should be noted that, due to the principle of special relativity, the Doppler equations, equations derived therefrom, and the machines, systems, methods and other techniques described in this application can be adapted to any embodiment or environment where an emitter and receiver are in motion relative to one another, for example: where the emitter is stationary and the receiver is moving, or where the receiver is moving and the emitter is stationary, or where the emitter and receiver are both in motion relative to a third frame of reference (e.g. the surface of the Earth, or a satellite in space). It shall be understood by one of skill in the art that the inventions disclosed herein are applicable to any frame of reference and to obtaining either the position of the receiver or the emitter, and that to merely change the frame of reference and/or the object for which the position information is sought does not take an embodiment outside the scope of the invention.

It also will be apparent to one of skill in the art that the background material, examples, descriptions, and other disclosures contained in this application are equally applicable to any inertial frame of reference and to finding the position of either the receiver or the emitter. Thus, for example, where this application discloses a method involving a stationary emitter and a moving receiver, one of ordinary skill in the art could, with routine algebraic manipulation, apply the method to a moving emitter and a stationary receiver. Likewise, where this application discloses a method involving calculating the position of the emitter relative to the position of the receiver, one of ordinary skill in the art could apply the method to calculate the position of the receiver with respect to the emitter.

Doppler-Based Derivatives

In certain embodiments, multiple independent measurements of Doppler-based derivatives may be made within the same short time interval typically used to obtain a single frequency estimate and these derivative-derived measurements may provide unique information concerning the position of the emitter. The instantaneous measurements are estimates of the observed Doppler shifted frequency (FOA), the rate at which the Doppler shifted frequency is changing (“FOA-DOT”), and Doppler curvature (“FOA-DOT-DOT”). The FOA equation is captured in the standard Doppler equation (Equation 1). The FOA-DOT measurement equation is the first derivative of the FOA equation with respect to time, and the FOA-DOT-DOT measurement equation is the second derivative of the FOA equation with respect to time. It has been found that these measurement equations (FOA, FOA-DOT, and FOA-DOT-DOT) may be used to predict the observed estimates based on hypothesized location and emitter frequency. The measurements can be obtained using any number of estimation techniques, including but not limited to polynomial fits to measured frequencies during a measurement window. Moreover, each of these measurements may provide unique information about the Doppler signature which can be used to constrain the emitter position. For example, FOA-DOT is much less sensitive to carrier frequency and also carrier drift than FOA.

It has also been found that measurements obtained from measurement windows separated in time can be combined to form new Doppler based measurements, which add new information to further increase the accuracy of the estimate of emitter location. For example, an FOA estimate may be obtained at time T1 and again at time T2. Taking the difference of these estimates results in a new measurement equation referred to herein as Difference of FOA (“DFOA”). The DFOA is different from FDOA in that it uses the observed Doppler-shifted frequencies collected on board a single moving platform as opposed to requiring simultaneous collections on board two separate platforms in order to estimate the standard cross-ambiguity function (“CAF”). Further, the DFOA measurement differs from FOA-DOT because of the separate measurement windows (i.e., differing time scales) over which the Doppler rate of change is estimated and due to the curvature of the Doppler shift.

In certain embodiments of the invention, the discrete difference measurement DFOA (“D-DFOA”) can also be used to approximate the curvature of the Doppler signature over time. Although the D-DFOA approximates the FOA-DOT-DOT over short time intervals, it contains different information which may improve geolocation estimates and, like DFOA and FOA-DOT, may become even more distinct as the time between samples used in the differencing increases. For example, the increase in information is most pronounced when the measurement windows span the period containing the zero Doppler measurement or the point of closest approach. Another combination is the difference of FOA-DOT at time 1 and FOA-DOT at time 2, which essentially provides a discrete time estimate of the FOA curvature.

The choice of the measurements to use will depend on the operational constrained optimizations (“CONOPs”). For example, for platforms targeting signals with high frequencies emitted by fast moving emitters, the instantaneous FOA, FOA-DOT and FOA-DOT-DOT, in conjunction with information on the known characteristics of the emitter (e.g. satellite trajectory data and/or satellite carrier frequency information) can provide the necessary information to solve for the position of the receiver (x_r, y_r) and the receiver velocity (v). If the altitude (z) is also required, these measurements can be combined with DFOA to provide the full x-y-z location of the receiver. As opposed to triangulation methods based solely on Doppler FOA measurements, the invention does not require waiting for a sufficient baseline to be established in order to triangulate the emitter location.

In certain collection scenarios there may be poor estimation accuracy over one or two measurement windows. In these scenarios, it may be necessary to employ several measurement windows in order to adequately measure the FOA and DFOA to support the desired geolocation performance. For example, when a collection platform collects low frequency signals from a relatively slow moving emitter, the Doppler derivatives may be smaller and measurement noise may limit the estimation accuracy. In this case, the FOA and DFOA may be measured over a longer period of time to insure that there is a measurable difference. This time interval, even for the slowest platforms, will provide much faster convergence than waiting for the baseline required to triangulate the receiver location based on FOA alone.

Once the measurements are collected, there are a number of methods that can solve these nonlinear equations for the state estimates. Note that some care must be taken not to fall into local minima when the number of measurements does not fully resolve the state ambiguities.

Referring once again to FIG. 1: this figure can also be used to illustrate a separate embodiment of the invention making use of the aforementioned Doppler based derivatives. A receiver 100 listens for high frequency signals transmitted by satellite 110 or 112. Due to the Doppler effect, when the satellite is approaching the receiver 110 the frequency of the signals observed by the receiver 100 will be increased. As the satellite approaches closer to the receiver 100, the difference between the source frequency and observed frequency will change as a function of the position of the satellite. Likewise, when the satellite is receding from the receiver 112 the frequency of the observed signals will be decreased; and as the satellite recedes further from the receiver, the difference between the source and observed frequency will continue to change as a function of the position of the satellite.

Referring to FIG. 4, a system according to one embodiment of the invention is illustrated. As shown, the system may comprise a receiver 460, which may be combined with or separate from an analog to digital converter. The system further comprises an FOA measurement system 420, an FOA-DOT measurement system 430, and FOA-DOT-DOT measurement system 440. Additionally, a processor or general computer 450 may be provided to calculate a non-linear least squares fit. As discussed above with respect to FIG. 1, each of the aforementioned systems may be on-board a moving or stationary platform. Likewise, the emitter system may be a satellite in a non-geosynchronous orbit. The individual measurements may be combined to determine a total solution, corresponding to the position of the emitter and having a particular confidence level.

Example 1

The performance of the above described DPF algorithm was examined for Doppler data simulated from a recorded UAV path, where the UAV attempts to determine the position of a 500 MHz RF emitter device in a 1 km² vicinity. Note that Example 1 illustrates simulated experimental data tracking an emitter with a platform having a known position and velocity. However, given that this experimental data is closely related mathematically to the problem of self-localization, the illustration provided by Example 1 is helpful in understanding self-localization. The UAV's flight path 501 over a collection period of 3 minutes (180 s) is shown in FIG. 5. As shown, the UAV starts at approximately [0;−600] x-y point with a height of 150 meters and travels at a relatively slow, constant velocity range of from about 8 to about 18 m/s (i.e., about 20-40 mph). The star 502 denotes the position of the target RF emitter, which is placed outside and about 0.5 km from the center of the UAV flight path.

Referring to FIG. 6, the simulated emitter carrier drift is illustrated 601. As shown, a drift of greater than about 20 Hz over the course of three minutes is well within the operational tolerances (40 parts per billion) for most commercial emitter devices.

Referring to FIG. 7, a comparison between the expected Doppler from a stable 500 MHz emitted signal 701 and the observed Doppler 702 from the emitter with the drift described in FIG. 6 and geometry of FIG. 5 is shown. The Doppler profile of the stable 500 MHz carrier 701 is shown for UAV flight path shown in FIG. 5. It will be appreciated that the Doppler shift for these speeds and geometry ranges from about −20 to about 20 Hz.

The Doppler profile of the drifting carrier 702 are also shown, minus the 500 MHz frequency reference. As shown, the drift changes the observed frequencies, resulting in a less clean signature of where the target emitter is located over the course of the collection. By tracking the emitter drift, a significantly more accurate and reliable geolocation can be obtained.

Referring to FIG. 8, the geolocation accuracy of 10 test runs of the DPF algorithm are compared to the FOA Cramer-Rao lower bound (CRLB), assuming an unknown but stationary carrier. As shown, each cross 801 represents one Monte Carlo solution of the particle solver evaluated every 20 seconds over a three minute flight path time period. The bold crosses 802 represent the mean of the 10 runs.

The results indicate that exemplary DPF algorithms may consistently obtain 50 meters or less geolocation errors after the first minute of flight and may be capable of closely approaching the FOA CRLB of about 20 m at the end of the three minutes. The initial input parameters to the particle filter solution were 11×11=121 particles spread uniformly over a 1000 m×1000 m grid and offset by 23 meters in each direction. The measurement noise used to simulate the FOA observables was 1 Hz². In certain embodiments, the DPF algorithm may use this sensor noise value for its Kalman filter's measurement covariance R. The Kalman's process covariance, Q, or the prediction variance for the frequency track, was 0.5 Hz and 0.05 Hz for the emitter carrier frequency and its rate of change, respectively. The line 803 denotes the FOA iterative CRLB. It was derived based on the assumption of zero process noise so that the bound can be efficiently found using the EKF provided the correct state information and in this case also the correct track of the emitter's carrier drift.

Referring to FIG. 9, the DPF algorithm and the Gauss Newton geolocation error is shown on a log 10 scale. The mean DPF algorithm error is shown as rectangular boxes 901, the FOA CRLB as circles 902, and the individual Gauss Newton trials as crosses 903. The Gauss Newton solution 903 assumes an unknown but stable carrier frequency. The logarithmic scaling shows the tendency of the Gauss Newton solution 1203 to vary widely and get trapped in local minima often as far as 10 km away from the actual solution and far outside the region of interest.

Example 2

Referring to FIG. 10, the performance of a DPF algorithm according to the invention is examined for the case where a target emitter 1302 is located inside the UAV's flight loop 1001. Note that Example 2 illustrates simulated experimental data tracking an emitter with a platform having a known position and velocity. However, given that this experimental data is closely related mathematically to the problem of self-localization, the illustration provided by Example 2 is helpful in understanding self-localization. As shown, the star 1002 denotes the position of the target RF emitter and the line 1001 shows the UAV's flight path over a collection period of 3 minutes (180 s). Similar to Example 1, the emitter comprises a variable 500 MHz RF carrier device having the same drift shown in FIG. 6. The UAV starts at approximately [0;−600] x-y point with a height of 150 meters and travels at a relatively slow velocity range of from about 8 to about 18 m/s (i.e., about 20-40 mph).

Referring to FIG. 11, a comparison between the expected Doppler from a stable 500 MHz signal 1101 and the observed Doppler 1102 from the emitter with the drift described in FIG. 6 and geometry of FIG. 10 is shown. It will be appreciated that the Doppler shift for these speeds and geometry ranges from about −20 to about 20 Hz. As shown, the drift changes the observed frequencies, resulting in a less clean signature of where the target is located over the course of the collection.

Referring to FIG. 12, the performance of the inventive DPF solution 1201 to the FOA CRLB 1202 given the new target location is shown. Again, the DPF algorithm error 1201 closely approaches the FOA CRLB 1202 with a mean value approaching 10 m accuracy for this target/sensor geometry. In fact, it will be appreciated that the preferred placement of the target inside the flight path dropped the expected error performance from about 30 meters to less than about 10 meters.

The Gauss Newton error 1203 is shown for individual trials. As shown, the Gauss Newton solution 1203 becomes trapped in local minima very far from the actual solution, while the DPF solution does not. These results demonstrate the capability of the particle filter approach to provide a robust, accurate FOA geolocation solution even for cases with substantial nonlinear emitter carrier drift.

Example 3

Finally, FIG. 13 compares the best potential performance of the two-sensor FDOA method 1301 to the single-sensor FOA method 1302 by comparing their respective CRLBs. Note that Example 3 illustrates simulated experimental data tracking an emitter with a platform having a known position and velocity. However, given that this experimental data is closely related mathematically to the problem of self-localization, the illustration provided by Example 3 is helpful in understanding self-localization. As shown, the geometry was the same as described in FIG. 5, except that the FDOA required a second collector placed at the location of [600;−600; 0]. Both bounds assumed a measurement error variance of 1 Hz². The results show that FDOA converges slightly faster, but that the performances of the two techniques at the end of the three minute collection are roughly comparable. Note also that the bounds assume that both FDOA sensors have absolute frequency measurements precise to within a Hz.

Unless specifically stated otherwise as apparent from the above discussion, it is appreciated that throughout the description, discussions utilizing terms such as “processing” or “computing” or “calculating” or “determining” or “displaying” or the like can refer to the action and processes of a data processing system, or similar electronic device, that manipulates and transforms data represented as physical (electronic) quantities within the system's registers and memories into other data similarly represented as physical quantities within the system's memories or registers or other such information storage, transmission or display devices.

The exemplary embodiments can relate to an apparatus for performing one or more of the functions described herein. This apparatus may be specially constructed for the required purposes, or it may comprise a general purpose computer selectively activated or reconfigured by a computer program stored in the computer. Such a computer program may be stored in a machine (e.g. computer) readable storage medium, such as, but not limited to, any type of disk including optical disks, CD-ROMs and magnetic-optical disks, read only memories (ROMs), random access memories (RAMs) erasable programmable ROMs (EPROMs), electrically erasable programmable ROMs (EEPROMs), magnetic or optical cards, or any type of media suitable for storing electronic instructions, and each coupled to a bus.

The exemplary embodiments described herein may be described as software executed on at least one computer or server, though it is understood that embodiments can be configured in other ways and retain functionality. The embodiments can be implemented on known devices such as a personal computer, a special purpose computer, cellular telephone, personal digital assistant (“PDA”), a digital tablet, an electronic gaming system, a programmed microprocessor or microcontroller and peripheral integrated circuit element(s), and ASIC or other integrated circuit, a digital signal processor, a hard-wired electronic or logic circuit such as a discrete element circuit, a programmable logic device such as a PLD, PLA, FPGA, PAL, or the like. In general, any device capable of implementing the processes described herein can be used to implement the systems and techniques according to this invention.

It is to be appreciated that the various components of the technology can be located at distant portions of a distributed network and/or the Internet, or within a dedicated secure, unsecured and/or encrypted system. Thus, it should be appreciated that the components of the system can be combined into one or more devices or co-located on a particular node of a distributed network, such as a telecommunications network. As will be appreciated from the description, and for reasons of computational efficiency, the components of the system can be arranged at any location within a distributed network without affecting the operation of the system. Moreover, the components could be embedded in a dedicated machine.

Furthermore, it should be appreciated that the various links connecting the elements can be wired or wireless links, or any combination thereof, or any other known or later developed element(s) that is capable of supplying and/or communicating data to and from the connected elements. The term module as used herein can refer to any known or later developed hardware, software, firmware, or combination thereof that is capable of performing the functionality associated with that element. The terms determine, calculate and compute, and variations thereof, as used herein are used interchangeably and include any type of methodology, process, mathematical operation or technique.

The invention described and claimed herein is not to be limited in scope by the specific embodiments herein disclosed since these embodiments are intended as illustrations of several aspects of the invention. Any equivalent embodiments are intended to be within the scope of this invention. Indeed, various modifications of the invention in addition to those shown and described herein will become apparent to those skilled in the art from the foregoing description. Such modifications are also intended to fall within the scope of the appended claims. All publications cited herein are incorporated by reference in their entirety.

Claims

1. A method of determining a location of a radio frequency (rf) receiver using an rf emitter located on a single moving platform, the method comprising:

A. receiving, by the receiver, an rf signal transmitted by the emitter over a period of time;

B. measuring, by the receiver, frequency of arrival (FOA) data of the rf signal over the period of time, the FOA data comprising a doppler-shifted frequency of the emitter, a location of the emitter, and a velocity of the emitter;

C. determining an estimate of a frequency measurement bias of the receiver based on the measured FOA data and an estimate of the location of the receiver, wherein the estimated frequency measurement bias of the receiver includes an estimate of a frequency drift; and

D. determining an updated estimate of the location of the receiver based on the estimate of the frequency measurement bias of the receiver.

2. A method according to claim 1 further comprising measuring initial FOA data over an initial time period and calculating the estimate of the location of the receiver by modeling a region of interest with a plurality of particles.

3. A method according to claim 2, wherein each particle represents a possible location of the receiver and is associated with a frequency drift conditioned on the location and the initial FOA data.

4. A method according to claim 2, wherein each particle represents a possible frequency measurement bias estimate and is associated with a location of the receiver conditioned on the frequency measurement bias estimate, initial FOA data, and the observed signal frequency of the beacon.

5. A method according to claim 2, wherein the estimate of the location of the receiver is based on a priori information or spread uniformly out in a grid-like fashion.

6. A method according to claim 1, wherein the estimate of the frequency measurement bias of the receiver is calculated with a particle filter.

7. A method according to claim 6, wherein the particle filter weights particles having a frequency drift that best corresponds to a model of an expected receiver frequency measurement bias drift.

8. A method according to claim 7, wherein the model is a process model of a localized extended Kalman filter.

9. A method according to claim 7, wherein the model is a parametric fit to the estimated frequency drift.

10. A method according to claim 7, further comprising resampling the weighted particles.

11. A method according to claim 10, wherein the updated estimate of the location of the receiver is determined based on a covariance of the resampled particles.

12. A method according to claim 1 comprising:

receiving, by the receiver, a second rf signal transmitted by the emitter over a second period of time;

repeating step C with the updated estimate of particle location to calculate an updated estimate of a frequency; and

repeating step D with the updated estimate of the frequency measurement bias of the receiver.

13. The method of claim 1, wherein the emitter is an orbiting satellite or aerial vehicle.

14. The method of claim 1, wherein the one or more of the group consisting of FOA-DOT, FOA-DOT-DOT, FDOA, and D-FDOA are used to calculate the updated estimate of the location of the receiver.

15. A system for determining a location of a radio frequency (rf) receiver comprising:

a receiver located on a platform, the receiver adapted to receive an rf signal transmitted by an emitter over a period of time; and

a processor in communication with the receiver adapted to:

measure frequency of arrival (FOA) data of the received rf signal over the period of time, the FOA data comprising a doppler-shifted frequency of the emitter, a location of the emitter, and a velocity of the emitter;

determine an estimate of a frequency measurement bias of the receiver based on the measured FOA data and an estimate of the location of the receiver, wherein the estimated frequency measurement bias of the receiver includes an estimate of a frequency drift; and

determine an updated estimate of the location of the receiver based on the estimate of the frequency measurement bias of the receiver.

16. A system according to claim 15 wherein the processor is further adapted to measure initial FOA data over an initial time period and calculate the estimate of the location of the receiver by modeling a region of interest with a plurality of particles.

17. A system according to claim 16, wherein each particle represents a possible location of the receiver and is associated with a frequency drift conditioned on the location and the initial FOA data.

18. A system according to claim 16, wherein each particle represents a possible frequency measurement bias estimate and is associated with a location of the receiver conditioned on the frequency measurement bias estimate and initial FOA data.

19. A system according to claim 16, wherein the estimate of the location of the receiver is based on a priori information or spread uniformly out in a grid-like fashion.

20. A system according to claim 15, wherein the estimate of the frequency drift is calculated with a particle filter.

21. A system according to claim 20, wherein the particle filter weights particles having a frequency measurement bias drift that best corresponds to a model of an expected frequency measurement bias drift.

22. A system according to claim 21, wherein the model is a process model of a localized extended Kalman filter.

23. A system according to claim 21, wherein the model is a parametric fit to the estimated frequency drift.

24. A system according to claim 21, wherein the processor is further adapted to resample the weighted particles.

25. A system according to claim 24, wherein the updated estimate of the location of the receiver is determined based on a covariance of the resampled particles.

26. The system of claim 15, wherein one of the group consisting of the emitter or the platform is an aerial vehicle or orbiting satellite.

27. The system of claim 15, wherein the processor employs one or more of the group consisting of FOA-DOT, FOA-DOT-DOT, FDOA, and D-FDOA to calculate the updated estimate of the location of the receiver.

28. A system according to claim 15 further comprising the rf emitter.

29. A system according to claim 15, wherein the processor is located on-board the platform.

30. A system according to claim 15, wherein the processor is located external to the platform and is in communication with the receiver via a network.