Multi-Target Tracking In The Littoral are not obvious at first glance.

It is evident from approaches

often adopted in current operational systems that these
Environment lessons (even some older ones) have not been absorbed by
the industry. Consequently this manuscript begins with a
Eric R. Bechhoefer, BFGoodrich brief review involving a very early fundamental choice of
Aerospace,Vergennes VT USA parametrization - and illuminates that choice in a way not
previously highlighted. The resulting implications provide a
James L. Farrell, Ph.D. basis for an expansion in scope that produced a new
Vigil Inc., Severna Park MD USA formulation for sea surface targets. Following that
explanation the formulation is presented, after which rigorous
simulation results are shown.
Abstract
Background
Existing state estimation algorithms for sea surface In the early 1970s an effort was begun, with the purpose
targets allow motion in three spatial directions. The approach of furthering the exploitation of computer capability in
described herein exploits directly all available information, airborne radars. Almost immediately in that effort, it was
with surface ships constrained to the geoid - producing a recognized (by one of the coauthors of this manuscript) that
model as close to reality as any representation humanly all earlier range track and angle track methods from previous
achievable. Not only is Earth curvature taken into account radars could be disregarded as obsolete. In one camp,
from the outset, but Earth oblateness and the geoid/ellipsoid Cartesian vector tracking was adopted with never a backward
height separation (taken into account by lookup table) are glance. Documentation of the approach [Refs. 1-6] was
sewed into the estimator algorithm, in the same form used for unfortunately obscure; this paper is intended to provide the
GPS data. These effects are included everywhere applicable, dissemination warranted.
not just as an external correction or an afterthought. To Concurrent with the earliest of events just described, an
enhance computational reliability still further, an alternate approach using spherical coordinates in a
appropriately stable (e.g., factorized) algorithmic approach is line-of-sight (LOS) coordinate frame was published [Ref. 7].
unhesitatingly adopted. Factorization provides another side That created a need to communicate some justification for
benefit of efficiency which, when added to the savings vector approach selection, e.g. the advantages of vector
inherent with fewer states,1 produces an almost irresistible tracking:
incentive to pursue this opportunity to reap full benefit from 1. High degree of linearity BOTH dynamics and
tracker input data received observables [Ref. 5]:
a. Correction of a common misconception: Vector
track observables can be expressed in a very nearly
Introduction
linear formulation.
2. Simplicity of implementation, both hardware and
The underlying commonality connecting several
software:
different tracking operations sometimes needs tailoring for a
a. Trackers from Refs [1-6] interface directly with INS
specific application. There are examples wherein one Kalman
output velocity only.
estimator formulation applies to multiple modes (e.g.,
b. Some trackers documented elsewhere require
air-to-air and air-to-ground with targets at short range), but
interface with gyros and accelerometers to achieve
there are counterexamples calling for some modification (e.g.,
comparable performance.
air-to-surface with ships at long range -- which is the primary
3. Isolation from stabilization error:
subject of this manuscript). Some of the reasons for these
a. Imperfect sensor pointing does not degrade
decisions go back several years (in fact decades), and many
estimation.
4. Versatility; applicability to multiple modes:
1 a. Adds, at most, computations involving outputs
Because the dominant matrix computations carry a cost only.
varying with the square or the cube of the dimensionality, 5. Extension to multitarget operation is straightforward.
savings from a smaller number of states can be more 6. Handoff to other sensor tracker is straightforward.
significant than initially expected. This issue becomes These items were frequently offered as justification for
important when the number of tracked objects grows large choosing a vector tracking approach, but one scenario
which quite commonly characterizes the sea-surface seemed to offer particularly dramatic persuasion: An
application. interceptor flying Eastbound at 1000 ft/sec tracks a target
flying Westbound at 1000 ft/sec on an antiparallel course where the angles are subtended. In both the Cartesian and
with 1000 ft separation distance. When the two aircraft reach the spherical coordinate representation initiall y cited, the
a sidelong position, the centripetal acceleration is (2000) * interceptor antenna phase center is the origin of the
(2000) / (1000) ft/sec/sec, over 120 g ! A Cartesian vector coordinate frame - this is an accepted convention. For
tracker will easily deduce a quiescent target velocity vector latitude and longitude a point near the geocenter {on the
but, with spherical coordinates, trackers will not be designed ellipsoid's evolute, to be precise; page 88 of Ref. [10]} is the
to follow accelerations anywhere near that level. Efforts to position reference. None of the advantages previously
follow extreme dynamics seemed especially futile when described are sacrificed by adopting geodetic latitude and
considering that range-channel acceleration is a fictitious longitude to express ship position - and some very
artifact, necessitated by rotation of the LOS coordinate frame. significant benefits are offered.
As if that plus the preceding itemization were not enough , an Some of the advantages provided by a geodetic
independent evaluation [Ref. 8] gave the Cartesian approach formulation are as immediate as they are obvious.
an edge in performance. No further consideration was given Compliance with the best Earth models available will be
to the approach in Ref. [7] -- despite its correctness from a automatically ensured by the following procedures:
theoretical standpoint. Reference [2] illustrated how, given 1. The state extrapolation part of the Kalman filter will use
Cartesian states and interceptor (ownship) velocity, estimates WGS-84
of familiar quantities (range rate, LOS rate, aspect angle, etc.) 2. Refinement with a data base for geoid height enables an
are easily computed on the basis of Kalman filter outputs. It intrinsic mean sea level (MSL) constraint everywhere; no
is now instructive to focus attention on formulation of target altitude states. The full span for target degrees of
dynamics in this manner. freedom in position and velocity can thus be covered by
First note that total target acceleration components are four rather than six states.
identical to the last three state variables - and thus need no The latter benefit can be significant - not because a
further processing whatever - if the acceleration "time Kalman filter imposes high computational demand {the
constant" from Ref. [9] is allowed to reach infinity {as in Eq. industry has come a long way since fixed-point computation
(5-64) rather than (5-63) of Ref. [10]}. Even with finite time necessitated dynamically scaled parameter values as in Ref.
constant, only a final-value adjustment is needed to scale the [2]} -- but because there can be a very large number of
target acceleration. That adjustment factor is based on gains tracked objects. To use the cube rule-of-thumb, the savings
used for velocity and acceleration states; much of the will be (4/6)3, a reduction by better than a factor of three.
industry appears unaware of (1) that factor's existence, (2) the Most important, performance is enhanced by enforcing geoid
need to apply it separately in along-range and cross-range model compliance -- not as an afterthought or a subsequent
channels, (3) superiority of performance achievable with an external adjustment, but as an inherent property of the
infinite time constant in many applications, (4) the need to formalism -- and everything humanly possible is being done
include a slow rotation of the target acceleration vector to force conformance of the model to the real world.
during extrapolation within the estimator algorithm (i.e., that
vector should be characterized as constant, not in an inertial Scenario
frame, but in a frame with one axis aligned to instantaneous Practical scenarios for applications under scrutiny can
total target velocity), and (5) acceptability of omitting effects be typified by 15000-ft radar altitude, creating a
of the adjustment just described on the state covariance range-to-horizon conforming to geometric mean {2 • altitude
matrix (though not the state vector) extrapolation {see page • Earth radius}1/2 and therefore an area of 2 π • altitude • Earth
180 of Ref. [10]}. radius = 50,000 nmi 2. To cite a rough figure a typical range /
With the Cartesian approach, formation of total target azimuth cell with 1 / 20 aspect ratio, corresponding to 300-ft
velocity is also quite straightforward, calling for mere vector and 1-nmi resolution, would imply a total of one million cells.
addition of the relative velocity to ownship velocity. All If one cell out of a thousand produced an above-threshold
target dynamics, then, can readily be formed from Kalman response, there could be up to 1000 active cells each time the
estimator outputs with minimal complexity imposed by area is scanned (i.e., 1000 • the fraction of cells not masked
ownship motion. It is actually this feature, rather than by landmass). A closer figure is obtained by omitting some
coordinate parametrization per se, that is key to the approach. area near the radar (which still retains most of the region to
Only when that principle is clearly understood can the be scanned) and dividing the remaining area into sectors ) for
significance of our next innovation be appreciated: a 2-deg beamwidth and a 360 deg/sec scan rate, there could
air-to-surface tracking of ships at long range should use be somewhat less than 180 • 20 • 130 (still nearly a half
latitude and longitude as position states. Immediately the million) cells. In that case a 0.001 detection ratio would
question arises: Isn't that backtracking or contradicting the produce up to almost 500 reports / sec (e.g., 400 with less
benefits just discussed at length? Answer: no - look at than 20% land masking).
 Patrol aircraft, such as the P-3 Orion, can conduct an Eqs. (1-28) and (7-31) in [Ref. 10] while holding the
open ocean search at altitude of 10,000 to 15,000 ft - giving a [unobservable] altitude estimate continuously at zero}, the
radar horizon of approximately 120 to 150 nautical miles. radar-to-target range vector is calculated
Reliable detection (e.g. Probability of detection) for a medium RT / O = RT - RO (1)
to small contact is significantly less, or around 105 nautical followed by the predicted measurement and H. The
miles under notional conditions (e.g. sea state 3, probability transformation To is calculated form (1-28) in Ref. [10] and:
of false alarm of 10-4, 1000 square meter vessel, APS -137 the range vector in observer's local coordinates is
radar in search mode). The test simulation used an initial R = To T RT / O , (2)
position of the target approximately 045 relative bearing at R = || R || , (3)
100 nautical miles. where ( ) T denotes the transpose operation.
In order to calculate the difference between the predicted
Mechanization sightline (e.g. bearing) and measured sightline to the target,
Track state vector consist of Latitude, Longitude, the radar sightline must be expressed in the same coordinate
Velocity North (in nautical miles per hour) and Velocity East system as R. This requires a coordinated scheme of track
(in nautical miles per hour), while the actual track update is stabilization, with identification of pertinent sightline
conducted in local Cartesian coordinates of the target (this directions from Ref. [5]:
allows multiple observers to exploit each others’
T True sightline running from the sensor to the
measurements, sharing data from multiple platforms).
target position
Coordinate frames defined for both geodetic position
and navigation references consist of three mutually I Indicated sightline as observed from an imperfect
orthogonal vectors chosen according to right hand sensor.
convention (I = J × K, J = K × I , K = I × J , I · I = J · J = K ·
K = 1 ). The frame rigidly attached to the earth is denoted as B The orientation of the sensor axis while the
(IE , J E , KE ). Conversion from geodetic position to local indicated sightline is being observed
coordinates is done via transformation into (IE , J E , KE ) using
(1-28) or (7-31) of Ref. 10, as appropriate. Specific notation is P Predicted sightline corresponding to an
as follows: anticipated direction based on a priori estimates
of the target state.
L Latitude
A Airframe reference.
l Longitude
G geographic reference.
VN North Velocity Familiar angular offsets are easily calculated from the
VE East Velocity differences between these unit vectors, e.g.,. measurement
error is I-T, sensor output is I-B, ideal sensor output T-B,
X target State Vector (L, l, V N , V E ) anticipated sensor output P-B, estimator input residual I-P.
The residual value is small, allowing for small angle
RT target position in Earth reference (IE , J E , KE ) approximation, essentially removing any nonlinearity in the
calculation of bearing (or elevation in a 3D system; for a more
RO observer position in Earth reference (IE , J E , KE )
complete description see Ref. 5.). With elevation information
∆t elapsed time from last measurement absent, there is a difference in the estimated bearing sightline
and the measured sightline, but this does not hinder
TO rotation Matrix to (IE , J E , KE ) from geographic projection (e.g., rotation) of the range vector RT / O into local
coordinates at observer's location coordinates. Bearing is the inverse tangent of the ratio of
East to North component of R (this is true for true and for
H measurement sensitivity matrix predicted values) , BRG = Arctan(R E / R N) and, for update,
For state extrapolatin between observations here, each the measurement matrix H is now efficiently calculated using
tracked object is maintained in its own target-centered a Jacobian of the range and bearing measurements with
reference (this improves both linearity and multisensor respect to range vector components;
operation, with all filtering performed in a local Cartesian
plane). Estimated position states are then:propagated
 RN RE 
according to the usual VN ∆t and sec (L) VE ∆t increments
 R R 0 0
H = R − RE (4)
for latitude and longitude respectively. Concurrently also N 0 0
maintaining Earth coordinates for radar and all targets {using  ( RN2 + RE2 ) ( RN + RE )
2 2

 An efficient mechanization of the Kalman filter, such as (1.5n 2+1.5n)m and a total number of multiplications of
Bierman {[Ref.11]} is used to filter in local coordinates giving (1.5n 2+5.5n)m and nm divisions where n is the number of
the best estimate for updated state. Using standard methods states, m is the number of measurements. This is compared
this is again transformed from geographic to (IE , J E , KE ) to to Kalman Stabilized Data Processing with the total number
maintain position concurrently in both Earth and geographic of additions of (4.5n 2+5.5n)m, total number of multiplications
conventions at all times. of (4.5n 2+7.5n)m and m divisios. Computation results show
this to be a good approximation.
Measures of Effectiveness APPENDIX A: TRACKING of AIRBORNE vs.
SEA SURFACE OBJECTS
Error is expressed here in terms of physical quantities
associated with the sensor. This will be through observation The following are not addressed within this immediate scope:
of the line of sight (LOS) vector, where the use of the angular 1. Passive (angles-only) tracking
displacement and range gives explicit insight into the track 2. Doppler observations from coherent radars
mechanization. Conceptually, we can characterize an angular 3. ECM, ECCM
offset error as the displacement of a dot relative to an 4. Trackable objects on land
intersection of cross hairs, referenced to the sightline from 5. Fusion with other sensors (e.g., FLIR, ESM, etc.)
observer to target . This is indicative of operation for any although supported, not necessary to tabulate here.
narrow field of view sensor (e.g. a few degrees), for which the 6. Centroiding
offsets are small angles. In essence these offsets are linearly 7. Landmass blanking
proportional to the azimuth (and, where applicable, elevation) 8. Association of each radar report with its own previous
residuals. Additionally, time derivatives of the range and reports
LOS azimuth are used to assess delay in the system. The last three items are quite relevant but, if covered
Range and bearing accuracies are determined by with any depth, would expand the scope too far for a single
calculations similar to the already presented track filter manuscript. Only a brief summary is attempted here.
mechanization expressions. In that mechanization the true First, radar reports are preprocessed by centroiding
range and bearing are of course unknown but, for evaluation echoes in adjacent cells (for this application, range and
of simulation results, T is known – thus the difference I-T is azimuth). Conceptually and qualitatively it makes sense to
used to calculate range and bearing error. Derivation of use a range-azimuth grid in this case but, at 100 miles, the
range rates is made by taking the dot product of the relative aspect ratio corresponding to 300 ft and 0.7 degree exceeds
velocity vector with the sightline direction. Bearing rate is 20:1. Thus a 5x5 grid will span twenty times as much distance
formed from the cross product of relative velocity and range in azimuth as it does in range; at 200 miles the figure of
vectors normalized by the range squared; difference obtained course changes to forty. It is widely recognized that, under
with true and with estimated values provides the error. these conditions, clustering behavior of tracked objects in
Computational effort, as stated before, can be decisive cells will vary dramatically with radar position. If a tracked
in applications with hundreds of tracked objects. object can straddle more than two azimuth cells for any
geometry, that same object can straddle well over five range
Results cells ) objects in many geometries can span far more range
The results of 60 Monte Carlo runs of 120 seconds show gates than azimuth cells. An increase in the number of range
only marginal improvement of the 4-state over the 6-state cells would reduce the imbalance in spatial resolution. As
filter. This could be the result of numerical issues. another ramification a "point" target could represent a large
Additionally, it is difficult to ensure that the plant noise for object whose long dimension is aligned within one AZ cell.
each filter is the same, as the 6-state filter effectively adds Finally it is noted that centroids are typically unweighted.
more power to the covariance matrix for a given spectral While it is realized that scintillation curbs the advantage of
density, which is used to estimate the unknown behavior of amplitude-weighting, it would not be unreasonable to allow
the target (see Figures 1 through 4). However, the real weights to depend to a limited extent on signal strength.
benefit is in computational efficiency; the 4-state has Radar return from landmass can be suppressed via a data
essentially a third the number of computations when base accurate to within radar resolution ) but uncompensated
compared to the 6-state. Additionally, when compared to coordinate errors exceeding resolution cell size could easily
Kalman Stabilized filter, the UDU factorization method is far overtax the processing load. Where accuracy cannot be
superior and numerically stabile (in UDU filtering, the guaranteed or where land/water boundaries introduce
covariance is maintained as an upper triangular format, where variations, conservatism is warranted. This raises the topic
positive definiteness is ensured). UDU Factorization of beam stabilization; some radars are without attitude data
methods {Ref. [11]} have a total number of additions of on the interface, in which case presence of azimuth rotation
can be detected by noting the time between consecutive Ohio, 1978.
North crossings of the radar beam. In any case, rotation of 4. Hedland and Farrell, "Simulation of Tracking Radar in
the radar-carrying aircraft must be compensated by one the Presence of Scintillation," NAECON Symposium,
means or another; absence of correction would not only Dayton Ohio, 1980.
degrade track accuracy but also interfere with landmass 5. Farrell and Quesinberry, "Track Mechanization
suppression. Alternatives," NAECON Symposium, Dayton Ohio,
If echoes from any tracked object are mixed with those 1981.
from any other object, the resulting time history will be 6. Farrell and Hedland, "Integrated Tracking Software for
fictitious. With large numbers of track files there clearly is Multimode Operation," AIAA Digital Avionics Systems
risk of confusion on a massive scale. Volumes have been Conference (DASC), Baltimore Md., Dec 1984.
written on this association problem and no attempt is made 7. Pearson and Stear, "Kalman Filter Applications in
here to duplicate all (or in fact any) of the methods invented Airborne Radar Tracking," IEEE TRANS Aerospace and
to alleviate the risk. This paper addresses the disposition of Electronic Systems, vol AES-10, May 1974.
radar data in each separate file, after correct associations 8. Asseo and Ardila, "Sensor-Independent Target State
have been made. Estimator Design and Evaluation," NAECON
Similarities and differences between the littoral and other Symposium, Dayton Ohio, 1982.
tracking environments are tabulated here after the plots.. 9. Singer, "Estimating Optimal Tracking Filter Performance
for Manned Maneuvering Targets," IEEE TRANS
APPENDIX B: REFRACTION EFFECTS on Aerosp & Electr Sys, vAES-6, July 1970.
10. Farrell, Integrated Aircraft Navigation, Academic
LONG RANGE TRACKS
Press1976 [now paperback only].
11. Bierman , Factorization Methods for Discrete Sequential
Fortunately, with the help of Reference [12] this issue can be
Estimation, Academic Press, 1978.
de-emphasized for cases wherein radar updates consist of
12. Bean and Dutton, Radio Meteorology .
only range and azimuth (not elevation) data. The reason is
traceable to the nature of the refraction phenomenon; Bean
and Dutton depict a geometric range, along an arc with
extremities at each end of the slant range segment, plus an
apparent radio range which accounts for speed variations
along that arc. To be succinct here, the top of page 360
characterizes the difference between slant range and apparent
radio range as "small compared to the height error" which
"constitutes over 95 percent of the total error" (bottom of
page 356; for a quick intuitive explanation consider the
difference in direction of the arc and the straight line slant
range). Range error can still be hundreds of feet; in modes
where this is excessive, a simple correction scheme should
suffice since range error is not the dominant effect of
refraction. No further effort was expended on this area but, Figure 1. Range Error, 4 vs. 6 State Filter.
before leaving the subject, it is worthwhile to recall the
benefit of our minimum-state estimation approach: ships are
constrained to the geoid throughout, ensuring derivation of
full benefit from the radar data to be received.

REFERENCES
1. Farrell and Quesinberry, Optimal Air-to-Air Tracking for
WX Radar (U)," Tri-Services Radar Symposium, West
Point N.Y., July 1974 (CONFIDENTIAL)
2. Farrell, Quesinberry, Morgan, and Tom, "Dynamic
Scaling for Air-to-Air Tracking," NAECON Symposium,
Dayton Ohio, 1975
3. Farrell, Tom, and Nemec, "Air-to-Air Designate/Track
with Time Sharing," NAECON Symposium, Dayton Figure 2. Range Rate Error, 4 vs. 6 State Filter
 Estimation approaches Optimal (Kalman) + suboptimal
(observers)

Number of states / axis Two or three (i.e., acceleration

optional)

State dynamics Linear

Observable / state relation Essentially linear if done

properly

Table 3: Differences, Littoral vs. Airborne Object Tracking

Characteristic Airborne Sea Surface
Figure 3. Bearing Error, 4 vs. 6 State Filter
Applicable clutter Moderate Very large
densities

3-D motion Unconstrained Constrained

Full instantaneous fix Rng / Az / El Rng / Az

content

Target speeds High Low

Target Several g Fractional g

maneuverability

Unintended Wind, gusting Waves,

disturbance source current

Req'd data rates > 3 Hz < 1 Hz

Figure 4. Bearing Rate Error, 4 vs. 6 State Filter (typical)

Typical allowable 1 sec. Tens of

Table 1. C o m p u t a t i o n a l Efficiency o f Various pre-track delay seconds
Implementations
Target State Full or partial Full
States Jacobian (FLOPS) Filter Update
Estimator (TSE)
(FLOPS)
coupling
4 State UDU 22 285
Influence of Earth Ext to TSE Ext + optional
6 State UDU 56 810
curvature Int
4 State Standard 22 913
algorithm Areas to be masked XMIT pulse overlap Land mass
6 State Standard 56 2685 + Main Beam
algorithm Clutter

Table 2. Similarities, Littoral and Airborne Object Tracking Customary origin Ownship Surface
location location

Applicable target densities Low to high; scenario-

dependent

Required accuracies Depends on usage

Achievable accuracies Dependent on sensor and

scenario

Search area Potentially wide

Association methods Same set (e.g., mult-hyp, min

accel, etc.)