Tracking method for a radar system

A tracking method for a signal echo system, including generating a plurality of gates for respective propagation modes on the basis of a target state prediction for a dwell time, and generating a target state estimate for the dwell time on the basis of target measurement points which fall within the gates.

FIELD OF THE INVENTION
 The present invention relates to a tracking method for a radar system, such
 as a phased array radar system or a bistatic radar system. Although the
 following discusses use for radar systems, the invention could also be
 applied to other signal echo systems, such as sonar systems.
 BACKGROUND OF THE INVENTION
 Radar signals returned from a target allow information to be determined
 concerning the slant range, azimuth and speed of a target relative to the
 receiving system of the radar system. The receiving system however
 normally receives a number of signals returned from the target which have
 different propagation paths or modes. Noise received by and induced in the
 receiving system can also be mistaken for a return signal from the target
 and needs to be taken into account. Tracking methods have been employed
 which track a target on the basis of signals relating to one propagation
 mode. Yet selecting one propagation mode neglects information relating to
 other modes which can be used to enhance the accuracy and sensitivity of
 the tracking method.
 BRIEF SUMMARY OF THE INVENTION
 In accordance with the present invention there is provided a tracking
 method for a signal echo system, including:
 generating a plurality of gates for respective propagation modes on the
 basis of a target state prediction for a dwell time; and
 generating a target state estimate for said dwell time on the basis of
 target measurement points which fall within said gates.
 The present invention provides a tracking method for a signal echo system,
 including:
 obtaining target measurement points for a dwell time;
 initiating tracking by obtaining an initial target state estimate from at
 least one of said points;
 determining a target state prediction for a subsequent dwell time on the
 basis of said target state estimate;
 generating a plurality of gates for respective propagation modes on the
 basis of the target state prediction; and
 generating a target state estimate for said subsequent dwell time on the
 basis of target measurement points for said subsequent dwell time which
 are within said gates.
 The target state estimate may be generated by applying association
 hypotheses to said measurement points in said gates and association
 probabilities to said hypotheses, obtaining conditional state estimates
 from the measurement points for each hypothesis and summing said
 conditional state estimates multiplied by said probabilities.
 The tracking initiating step can advantageously be performed for a
 plurality of propagation modes to initiate a plurality of tracking filters
 by generating a plurality of said target state estimates for said
 subsequent dwell time.
 The present invention further provides a tracking method for a signal echo
 system, including extending a target state vector to include additional
 parameters associated with a plurality of propagation modes, and
 accounting for measurement uncertainty associated with propagation path
 characteristics for said modes when updating target state estimates.

DETAILED DESCRIPTION OF THE INVENTION
 Bistatic radar systems employ separate transmitter and receiver sites, and
 include Over The Horizon Radar (OTHR) systems which direct transmission
 signals above the horizon for refraction by the ionosphere, known as
 skywave systems. OTHR systems also include surface wave radar systems
 which propagate radar waves along the surface of saltwater, and rely on
 the receiving system being able to detect objects by the radar signals
 reflected therefrom.
 An OTHR system 2, as shown in FIG. 1, includes a receiving system 4 and a
 transmitting system 6. The transmitting system 6 comprises an in-line
 array 7 of transmitting antennas located at the transmitter site and a
 control system 10 for feeding electrical signals to the antennas. The
 receiving system 4 comprises an in-line array 12 of receiving antennas and
 a control system 16 for processing the signals received by the antennas,
 which are located at the receiver site. OTHR systems include the Jindalee
 Facility in Alice Springs (JFAS) and the U.S. Navy's ROTHR system.
 The broad transmitting beam of the radar is directed towards areas of the
 ionosphere from which refracted signals are redirected to monitor a target
 3. The beam is effectively directed to a region or area in which a target
 is located. A number of targets may be located in one region and the
 receiver control system 16 is able to divide the energy returned from the
 illuminated region into a dozen smaller beams which can then each be
 divided into a plurality of range cells that are characterised by a
 respective distance from the receiving system 4. This allows the receiving
 system 4 to track a number of targets which are located in the illuminated
 region. The receive beams can also be divided into a plurality of velocity
 cells characterised by an object's velocity relative to the receiving
 system 4. This allows targets to be separated on the basis of their
 velocity if they cannot be separated on the basis of their distance from
 the receiving system 4. The transmitting and receiving beams can be moved
 or swept in synchronism, through a number of beam steer positions, with
 the time being spent at any given position being referred to as the dwell
 time. Measurements obtained from the radar signals or echoes received
 during each dwell time are referred to as dwells.
 The control software of the control system 16 is able to obtain four
 parameters pertaining to a target from each dwell, and these are the
 propagation path length or slant range (R), azimuth (A), Doppler frequency
 shift or radial speed (D) and signal strength based on a signal to noise
 ratio (SNR) measurement. These are referred to as the RAD or radar
 coordinates. The set of measurements from a dwell also includes clutter
 and detections from other targets.
 The dwells can be graphically represented by plotting them as candidate
 detection points on a three dimensional axis, as shown in FIG. 2, for
 dwell t=k, where one axis represents R, the other A and the third the D
 values. For any dwell time t=k of the order of 100 or 1000 candidate
 detection points 50 may be determined by the receiving system 4. Some of
 the points 50 may correspond to a target and others may simply relate to
 clutter echoes or noise intrinsic in the transmitting system 6 or
 receiving system 4. Clutter echoes arise from backscatter from the ground
 or objects which are not of interest, such as meteors. The OTHR system 2
 is also subject to multipath propagation in that there is more than one
 single path for echoes returned from a target due to a number of different
 ionospheric layers 54 at different heights 53 which refract echoes down to
 the receiving system 4. as shown in FIG. 3. There may be up to four
 different reflecting layers F.sub.0, F.sub.1, F.sub.2 and F.sub.3
 resulting in several echoes returned from a target, corresponding to
 reflections from combinations of these layers. Propagation modes are
 described by the layers from which the signal is refracted. For example,
 F.sub.0 -F.sub.1 is the propagation mode for a transmit path via layer
 F.sub.0 and a receive path via layer F.sub.1, where T represents the
 target 3, as shown in FIG. 3. Whilst the propagation path for a candidate
 detection point 50 is not known, the height of the different layers can be
 determined using commercial ionospheric sounders which provides some
 information concerning the relationship between points of different
 propagation modes for the same target. Knowing the heights and properties
 of each layer gives an indication as to expected RAD measurements of
 different propagation modes.
 The state of the target, at a given dwell k, can be represented by
 ##EQU1##
 where r is the range, a the azimuth, r the range rate and a the azimuth
 rate. Equations of motion can be used to describe the target dynamics, for
 example, a constant velocity target would, if the time T between dwells
 were constant, obey
EQU r(k)=r(0)+rkT
EQU a(k)=a(0)+akT. (2)
 This can be expressed in known state-space form as
EQU x(k+1)=F(k).times.(k)+v(k) (3)
 where F(k) is a known matrix, for instance in the case of a constant
 velocity target
 ##EQU2##
 where T.sub.k is the time between dwells k and k+1. The term .nu.(k)
 represents zero-mean, white Gaussian process noise as used in standard
 Kalman filtering. The covariance matrix Q(k) of .nu.(k) is assumed to be
 known.
 A currently used tracking method, based on the probabilistic data
 association (PDA) filter, as described in Y. Bar-Shalom and T. E.
 Fortmann, "Tracking and Data Association", Academic Press, 1988. performs
 tracking in the radar coordinates R, D, A, A, as illustrated in FIG. 2. A
 track is initiated by selection of a single, noisy measurement 50 with the
 unknown azimuth rate A being initially set to correspond with a
 hypothesised azimuth crossing rate, usually zero. Further measurement
 selection is accomplished by taking only those measurements which fall
 inside a validation gate 70 around the next expected position of the
 target measurement. This method does not require knowledge of the mapping
 of radar to ground coordinates during tracking. A disadvantage of this
 method is that it fails to use the information conveyed by multiple
 detections arising from multipath propagation. Also the presence of
 multipath propagation may cause multiple tracks 60, 62 and 64 to be
 generated for a single target, as shown in FIG. 4, when tracking is
 performed in radar coordinates using conventional filters such as PDA. If
 the tracks 62 and 64 closely conform with the expected separations for the
 hypothesised modes, they can be considered to relate to the same target 52
 of FIG. 5, whereas a track 66 which diverges excessively can be dismissed
 as corresponding to another target or to clutter. Such a situation,
 commonly arising in conventional PDA tracking, requires a fusion or
 clustering operation to group multimode tracks pertaining to the same
 target together. This allows a track to be identified with a particular
 propagation mode. A further stage of coordinate registration is then
 required to map the tracks to ground coordinates for geographical display
 to the radar operators.
 The preferred embodiment described herein uses explicit knowledge of the
 ionospheric structure including virtual heights, as provided by
 ionospheric sounders or by other means, to account for and take advantage
 of multipath propagation during track initiation and tracking. This is
 distinct from conventional approaches which only expect a single detection
 per target and are unable to benefit from the additional target-related
 information conveyed by multipath detections. The gain in tracking
 performance arising from multipath detections of a single target is
 important when the probability of target detection via some or all of the
 various propagation modes is low.
 The target state is taken to be as in equation (1), where r is the ground
 range 8 across the surface of the earth 9, a is the true azimuth, r is the
 ground range rate and a is the true azimuth rate. The true azimuth a is
 the complement of the angle a, which is the angle between the projected
 ground range r and the axis of the receiver array 12, as shown in FIG. 5,
 i.e. a=(90.degree.-a). Tracking is performed in ground coordinates,
 although other frames of reference, for instance a preferred propagation
 mode, may be used to describe the target dynamics and relate these to the
 other measurement coordinates.
 The conversion between the ground and radar coordinates can be represented
 as
 ##EQU3##
 where, at time k, R is the measured slant range, A the measured azimuth, D
 the Doppler speed (slant range rate), h.sub.r the virtual ionospheric
 height 53 on the receive path, and h.sub.t the virtual ionospheric height
 54 on the transmit path, as shown in FIG. 5. The slant range R may be
 defined as one half of the total path length from the transmitter 7 via
 the target 52 to the receiver 4. The measured azimuth or coning angle A is
 the complement of the angle A between the incoming ray 57 and the axis of
 receiver array axis 12. The Doppler speed D is proportional to the rate of
 change of the total path length.
 The various propagation modes can be labelled according to the
 corresponding outbound and return propagation mode combination F.sub.0
 -F.sub.0, F.sub.0 -F.sub.1, F.sub.1 -F.sub.0, . . . , F.sub.2 -F.sub.2 for
 a target 52, as shown in FIG. 3. For four possible ionospheric layers
 F.sub.0, F.sub.1, F.sub.2, F.sub.3 with heights h.sub.0, h.sub.1, h.sub.2,
 h.sub.3, these modes may be numbered from 1 to 16 respectively. Hence we
 may write the measurement process for the various propagation modes in
 terms of the target state x(k) as
 ##EQU4##
 where H.sub.1 (x(k))=H(r(k), a(k), r(k); h.sub.0, h.sub.0), H.sub.2
 (x(k))=H(r(k), a(k), r(k); h.sub.0, h.sub.1), etc., and the assumed number
 of possible propagation modes n may vary with time. In the above, w.sub.i
 (k) is a zero-mean, white Gaussian sequence with known covariance R.sub.i
 (k) representing the assumed measurement noise terms. The actual form of
 the non-linear measurement functions H.sub.i (.multidot.) above is
 determined by the geometry of the ionospheric model as shown in FIG. 5,
 and will depend on the virtual heights of the ionospheric layers h.sub.r
 and h.sub.t 53 and 55, and the location and separation of the receiver and
 transmitter arrays 7 and 12, among other factors.
 Since the virtual ionospheric heights h.sub.i in FIG. 5 may only be
 approximately known, but are assumed to vary slowly in comparison with the
 target dynamics, they can be included in the state vector x(k) and
 estimated along with the dynamical variables describing the target. In
 this case we have instead of equation (1)
 ##EQU5##
 with each virtual height satisfying an equation of the form
EQU h.sub.i (k+1)=h.sub.i (k)+.nu..sub.i (k) (8)
 where .nu..sub.i (k) is a small process noise term.
 Converting to the ground frame of reference requires the selection of an
 outbound and return propagation mode combination F.sub.t and F.sub.r with
 corresponding virtual heights h.sub.t and h.sub.r.
 The inverse transformation to equation (5) can be represented by
 ##EQU6##
 and follows from the assumed geometry indicated in FIG. 5.
 Hereinafter the state prediction and associated prediction covariance are
 denoted by x(k.vertline.k-1) and P(k.vertline.k-1) and an updated state
 estimate and state error covariance are denoted by x(k.vertline.k) and
 P(k.vertline.k).
 At some arbitrary time 0, tracking is initiated by selecting an initial
 point 50 which may correspond to a hitherto unobserved target. Since the
 propagation mode which gave rise to this measurement is a priori unknown,
 an initial target state estimate x(0.vertline.0) for equation (1) cannot
 be inferred from equation (9) unless a given propagation mode, or
 equivalently the ionospheric heights for the transmit and receive paths,
 is assumed. The preferred method is therefore to initialise n tracking
 filters, one for each possible initial propagation mode. Each filter
 assumes a particular initial propagation mode with corresponding virtual
 heights h.sub.r and h.sub.t in order to assign its initial state estimate
 using equation (9) based on the first measurement point 50. The estimate
 of the initial target azimuth rate is set to some starting value, usually
 zero. An initial state error covariance P(0.vertline.0) is also assigned
 and is taken to be large enough to cover the initial uncertainty in target
 position and velocity. Other methods of initialisation are possible using
 data from more than a single radar dwell; but the previously described
 method is the simplest among these. Of the n filters initiated from the
 measurement 50, the filter based on the correct initial propagation mode
 assumption can be expected to perform the best and thus its state
 estimates would be more accurate (in the sense of having smaller errors on
 average) than those of the other filters initiated with it. As the
 processing proceeds, it becomes clear by observation of the state
 estimates which, if any, of the n filters initiated as above is compatible
 with a target whose dynamical model is assumed to be as expressed in
 equation (2).
 The recursive processing required by each tracking filter, initialised as
 above, is now described. The aim of the processing is to compute, in a
 recursive manner, approximate conditional mean x(k.vertline.k) and
 covariance P(k.vertline.k) estimates of the target state, based on the
 measurement data, including virtual ionospheric height measurements, up to
 time k, Y(1), . . . , Y(k), where Y(i) represents the set of measurements
 received in dwell i. The estimated target track is provided by plotting
 the range and azimuth values from x(k.vertline.k). The accuracy of the
 track is indicated by the size of the standard deviations which can be
 obtained from the state error covariance P(k.vertline.k).
 The dynamical target model of equation (3) is used to predict where each
 measurement would appear during the next dwell in the absence of
 measurement noise under each propagation mode. The state prediction
 x(1.vertline.0) at time 1 is given, in the usual manner of Kalman
 filtering, as described in Y. Bar-Shalom and T. E. Fortmann, "Tracking and
 Data Association", Academic Press, 1988, as
EQU x(1.vertline.0)=F(0)x(0.vertline.0) (10)
 with associated covariance
EQU P(1.vertline.0)=F(0)P(0.vertline.0)F'(0)+Q(0) (11)
 where F' is the transpose of the transition matrix F in equation (3).
 Instead of generating one gate 70, x(1.vertline.0) is used to generate n
 gates 72, 74 and 76 in the measurement space for each tracking filter,
 corresponding to the respective propagation modes F.sub.0 -F.sub.0,
 F.sub.0 -F.sub.1, . . . , etc., as shown in FIG. 6. The measurement
 predictions for the respective propagation modes are therefore
 ##EQU7##
 The associated measurement prediction covariances are
 ##EQU8##
 where J.sub.i (1) is the Jacobian matrix of the non-linear measurement
 function H.sub.i (.multidot.) in equation (6) evaluated at the state
 predictions x(1.vertline.0). The validation gate for each propagation mode
 is an ellipsoidal region in RAD space defined by
EQU G.sub.i (1)={y.epsilon.IR.sup.3 :[y-y.sub.i (1.vertline.0)]'S.sub.i
 (1).sup.-1 [y-y.sub.i (1.vertline.0]&lt;.gamma..sub.i } (14)
 where .gamma..sub.i, defines the size of the validation gate. The
 probability that a target falls inside the gate i is denoted
 P.sub.G.sup.i, while the probability of detecting a target via the ith
 propagation mode is denoted P.sub.D.sup.i. This is illustrated in FIG. 6
 for three gates 72, 74 and 76 centred on measurement predictions 82, 84
 and 86 for three propagation modes F.sub.0 -F.sub.0, F.sub.0 -F.sub.1, . .
 . , etc. The gates may or may not overlap. The validation region is
 defined as the union of the validation gates or some region which includes
 their union. Points which fall inside the validation gates are accepted as
 possibly relating to the target 52 and are used together with the state
 prediction x(1.vertline.0) in order to update the state estimate
 x(0.vertline.0) to yield x(1.vertline.1). The corresponding state error
 covariance is also updated to P(1.vertline.1). This process is recursive
 and can be represented as follows:
EQU x(0.vertline.0).fwdarw.x(1.vertline.0).fwdarw.x(1.vertline.1).fwdarw.x(2.ve
 rtline.1).fwdarw. . . .
EQU P(0.vertline.0).fwdarw.P(1.vertline.0).fwdarw.P(1.vertline.1).fwdarw.P(2.ve
 rtline.1).fwdarw. . . . (15)
 The state estimate x(k.vertline.k) is an approximate
 minimum-mean-square-error estimate of the target state x(k) based on all
 the information 48 from dwells 0 through k of the form given in FIG. 2
 including multiple detections of the same target due to multipath
 propagation. The estimate is approximate because it assumes that the
 probability density function of the true target state is Gaussian
 conditioned on all the measurement data.
 To determine the updated target state x(k.vertline.k) and its covariance
 P(k.vertline.k), the measurements falling within the gates 72, 74, 76.
 etc. are used in a probabilistic data association framework as described
 in Y. Bar-Shalom and T. E. Fortmann, "Tracking and Data Association",
 Academic Press, 1988. which, in addition to consideration of a measurement
 being from a target or due to clutter, includes association hypotheses for
 the possible propagation modes which may have produced the measurements. A
 target existence or confidence model is also incorporated in the filter as
 described in S. B. Colegrove, A. W. Davis and J. K. Ayliffe, "Track
 Initiation and Nearest Neighbours Incorporated into Probabilistic Data
 Association", J. Elec. and Electronic Eng., Australia, Vol. 6. No. 3, pp.
 191-198, 1986, to aid in track maintenance (confirmation, deletion, etc.).
 The probability that the target exists at time k given data to time k is
 denoted P.sub.E (k.vertline.k). Target existence is modelled as a 2-state
 Markov Chain so that the predicted probability of target existence P.sub.E
 (k.vertline.k-1) satisfies
EQU P.sub.E (k.vertline.k-1)=.DELTA..sub.0 P.sub.E
 (k-1.vertline.k-1)+.DELTA..sub.1 {1-P.sub.E (k-1.vertline.k-1)} (16)
 where the two transition probabilities .DELTA..sub.0 and .DELTA..sub.1 are
 defined by
EQU .DELTA..sub.0 =Pr(target exists at time k.vertline.target exists at time
 k-1)
EQU .DELTA..sub.1 =Pr(target exists at time k.vertline.target does not exist at
 time k-1)
 An arbitrary initial value of P.sub.E (0.vertline.0)=0.5 is assumed.
 As an illustration of the filtering procedure, consider gates 72 and 74
 associated with propagation modes F.sub.0 -F.sub.0 and F.sub.0 -F.sub.1,
 and gate 76 associated with propagation mode F.sub.1 -F.sub.1, with
 respective centres given by the measurement predictions y.sub.1
 (k.vertline.k-1), y.sub.2 (k.vertline.k-1) and y.sub.3 (k.vertline.k-1),
 82, 84 and 86. We will number these propagation modes as 1, 2 and 3,
 respectively when referring to the measurement predictions. Suppose that
 the gate 72 contains two measurements y.sub.1, Y.sub.2 90 and that gate 74
 contains one measurement y.sub.3 92, while gate 76 does not contain any
 measurements. The 7 association hypothesis (numbered from -1 to 5) which
 can be applied are:
 (-1) The target does not exist.
 (0) The target exists but all validated measurements y.sub.1, y.sub.2 and
 y.sub.3 are clutter.
 (1) y.sub.3 and y.sub.3 are clutter, y.sub.2 is a target detection via
 propagation mode F.sub.0 -F.sub.0.
 (2) y.sub.2 and y.sub.3 are clutter, y.sub.1 is a target detection via
 propagation mode F.sub.0 -F.sub.0.
 (3) y.sub.3 is a target detection via F.sub.0 -F.sub.1, and both y.sub.1,
 and y.sub.2 are clutter.
 (4) y.sub.3 is a target detection via F.sub.0 -F.sub.1, y.sub.1 is a target
 detection via F.sub.0 -F.sub.0 and y.sub.2 is clutter.
 (5) y.sub.3 is a target detection via F.sub.0 -F.sub.1, Y.sub.2 is a target
 detection via F.sub.0 -F.sub.0 and y.sub.1 is clutter.
 For each of the possible associated hypothesis above, a conditional target
 state estimate x.sub.i (k.vertline.k) can be formed from the predicted
 state estimate x(k.vertline.k-1) using the extended Kalman filter theory,
 as described in G. W. Pulford and R. J. Evans, "Probabilistic Data
 Association for Systems with Multiple Simultaneous Measurements",
 Automatica, Vol. 32. No. 9. pp. 1311-1316, 1996. Omitting some time
 indexes and writing x=x(k.vertline.k-1) for equation (10),
 P=P(k.vertline.k-1) for equation (11), and y.sub.i =y.sub.i
 (k.vertline.k-1) for equation (12), the conditional state estimates in
 this case are given by
 ##EQU9##
 where the terms J.sub.i ; and S.sub.i are as in equation (13). The
 corresponding conditional state error covariances P.sub.i (k.vertline.k)
 are given by
 ##EQU10##
 where c.gtoreq.1 is a scaling factor reflecting the increased uncertainty
 in the case that the target does not exist.
 The computation of the probability of each associated hypothesis, called
 the association probability, can be illustrated by assuming uniformly
 distributed clutter measurements in the radar measurement space, and a
 Poisson model, as described in Y. Bar-Shalom and T. E. Fortmann, "Tracking
 and Data Association", Academic Press, 1988, with spatial density .lambda.
 for the number of clutter points inside the validation region. Also, we
 let the probability of target detection P.sub.D via any propagation mode
 be identical, and the gate probabilities P.sub.G be identical. P.sub.D and
 P.sub.G are parameters that are given values which are selected to extract
 optimum performance from the tracking method given the operating
 characteristics of the system 2. The total volume of the validation gates
 at time k is V.sub.k. If two or more of the gates overlap, V.sub.k can be
 approximated, for instance, as the volume of the largest gate. The
 association probabilities .beta..sub.i (k), defined as the probability of
 the respective association hypothesis i conditioned on all measurement
 data up to the current time k, can be expressed as described in G. W.
 Pulford and R. J. Evans, "Probabilistic Data Association for Systems with
 Multiple Simultaneous Measurements", Automatica, Vol. 32. No. 9. pp.
 1311-1316, 1996, by:
 ##EQU11##
 where N{y; y, S} is a multivariate Gaussian density iny with mean y and
 covariance S, and .delta.(k) is a normalisation constant, chosen to ensure
 that the association probabilities sum to unity.
 The updated target state estimate for the filter is obtained by summing the
 conditional state estimates with weightings determined by their respective
 association probabilities as
 ##EQU12##
 The state error covariance P(k.vertline.k) is obtained using standard
 techniques from Gaussian mixtures described in Y. Bar-Shalom and T. E.
 Fortmann, "Tracking and Data Association", Academic Press, 1988 as
 ##EQU13##
 The updated probability of target existence P.sub.E (k.vertline.k) is
 obtained as
EQU P.sub.E (k.vertline.k)=1-.beta..sub.-1 (k). (22)
 Track maintenance is achieved by thresholding the target existence
 probability according to
 ##EQU14##
 where P.sub.DEL and P.sub.CON are the track maintenance thresholds. Note
 P.sub.DEL &lt;P.sub.CON. Since P.sub.E (k.vertline.k) may vary considerably
 from dwell to dwell, it is better to use the average value of P.sub.E
 (k.vertline.k) over the last few dwells for track confirmation in equation
 (23).
 The above method is easily extended to arbitrary numbers of measurements
 falling inside the validation gates and to arbitrary numbers of
 propagation modes. Arbitrary clutter probability density functions and
 non-identical gate and detection probabilities are also easily
 accommodated within this framework.
 A tracking filter as described above has been implemented in software using
 the C programming language and executed on a Digital Equipment Corporation
 175 MHz Alpha workstation. The preferred implementation assumes 4
 propagation modes corresponding to F.sub.0 -F.sub.0, F.sub.0 -F.sub.1,
 F.sub.1 -F.sub.0 and F.sub.1 -F.sub.1. The virtual heights of the F.sub.0
 and F.sub.1. ionospheric layers are included as state variables in
 equation (7) and these are estimated from noisy measurements along with
 the range, azimuth, range rate and azimuth rate of the target.
 Many modifications will be apparent to those skilled in the art without
 departing from the scope of the present invention as herein described with
 reference to the accompanying drawings.