A State Space Technique for Wildlife Position Estimation Using NonSimultaneous Signal Strength Measurements
Abstract
\parttitleBackground Fixed veryhighfrequency (VHF) antenna arrays have been a standard technique for remote tracking of radiotagged wildlife since the 1960s. In recent years, a growing network of coordinated, fixed VHF arrays on a shared frequency has expanded throughout the Western Hemisphere, but our ability to estimate animal trajectories at regional scales is limited by the fact that detections consist of irregular time series of signal strength values as tagged individuals move in and out of detection range of antenna beams within the array. Therefore more advanced modeling methods are needed to refine estimates of animal locations. A novel statespace technique is presented to estimate the location and airborne movements of VHF tagged wildlife individuals with fixed VHF arrays. The approach combines a movement model (Ornstein Uhlenbeck random process in the transverse (horizontal) plane and a Cox IngersollRoss process in the vertical direction) to ensure biologicallyconsistent trajectories in threedimensions, and an observation model to account for the effect of range, altitude and bearing angle on the received signal strength. The observation model of received signals accounts for lowend saturation from receiver noise, highend saturation from receiver nonlinearities as well as a wireless multipath phenomena, which modulates the received signal according to range, the altitude and radiation characteristics of the Yagi array. A pattern function for the Yagi array is synthesized that facilitates linearization of the received signals and subsequent application of Kalman filtering.
\parttitleResults We first validate the model using a simulated trajectory and then estimate the spacetime trajectory of a migrating VHFtagged shorebird, which was tracked with a regional automated radio telemetry network. The algorithm accurately predicted the average movement trajectory given the system parameters and the initial conditions (average error 1 km). The modeled shorebird track represents a first estimate in threedimensional (3D) of a radio tagged bird using a fixed telemetry array, and was qualitatively reasonable, but exhibited some sensitivity in the vertical plane and to initial conditions.
\parttitleConclusions From an applied ecology perspective, we feel that our technique is flexible and broadly applicable to other telemetry networks and ecological systems. Thus, given appropriate array configuration and modification of the observation and process models, our technique could be applied to the analysis of automated VHF telemetry data across a wide range of spatial and temporal scales, from sitespecific studies using targeted arrays, to coordinated digital VHF tracking efforts that span the Hemisphere.
addressref=aff1, corref=aff1, email=janaswamy@ecs.umass.edu, pamela.h.loring@eco.umass.edu, JamesMcLaren@cunet.Carleton.ca ]\initsRJ\fnmRamakrishna \snmJanaswamy addressref=aff2, email=pamela.h.loring@eco.umass.edu ]\initsPL\fnmPamela \snmLoring addressref=aff3, email=JamesMcLaren@cunet.Carleton.ca ]\initsJM\fnmJames D. \snmMcLaren
Radio telemetry \kwdstatespace \kwdKalman filtering \kwdOrnstein Uhlenbeck process \kwdCoxIngersollRoss process \kwdbird migration \kwdradio multipath \kwdaeroecology
1 Background
Accurate estimation and assessment of animal movements is an essential tool towards understanding and conserving animal populations and ecosystem processes [1, 2, 3], [4]. Radio telemetry has been a standard technique for tracking wildlife since the 1960s [5], and to date remains one of the sole options for collecting individualbased tracking data of smallbodied ( 100 g) species [2]. For example, accurate positioning of birds fitted with lightweight ( g) VHF transmitters is possible using handheld receivers; manual tracking is however limited by effort, available transportation and knowledge of an animal’s whereabouts following release [5, 6, 7].
The recent expansion of coordinated digital VHF telemetry using automated radio telemetry systems has enabled a wide range of wildlife taxa to be tracked at unprecedented spatial and temporal scales. Automated radio telemetry systems can provide roundtheclock detection for relatively long periods (e.g. 4+ months for 1g units with a 5second pulse rate) whenever tagged individuals are within detection range of receiving antennas [8, 9]. Such systems have enabled studies to be conducted in a wide range of environments, from grassland [10] to tropical rainforest [11], to track a variety of taxa including reptiles [12, 13], birds [14, 15, 16, 17, 18], terrestrial mammals [19], and bats [20]. These studies typically base location on triangulation via implied bearings, estimated in turn via relative signal strengths between multiantenna array elements. More recently, more accurate methods have been developed using time of arrival and highprecision clocks [21], [22]. These studies have all been applied to ranges within 520 km and animals either on the ground or presumed to be flying at low altitudes. Contrastingly, a coordinated network of automated radio telemetry stations, the Motus Wildlife Tracking System, was piloted in northeastern North America in 2012 and has since expanded throughout western Hemisphere (www.motuswts.org). These stations have recorded over 250 million detections from more than 9,000 individuals representing 87 species of birds, bats, and insects tracked with digital VHF transmitters on a shared frequency [23]. However, given the data collected by automated radio telemetry stations typically consists of unevenlyspaced time series of signal strength values received by single beams or sporadically by multiple antenna beams, estimating the true positions and trajectories of radiotagged animals involving changes in flight altitude at regional scales remains a quantitatively complex and an unresolved area of research [24, 23].
2 Methods
2.1 Aim, Design and Setting
The objective of our study is to provide a framework to estimate the position of radiotagged birds based on timeseries data logged by an automated radio telemetry array, which accounts for the 3D structure. In this study, we used digitally coded VHF transmitters (NTQB42; Lotek Wireless; hereafter nanotag) that weighed 1 g and measured mm, with an external 18 cm long wire antenna. Each transmitter was programmed to emit signals at 166.38 MHz on a pulse duration of 5.3 s, from activation through the end of battery life (approximately 170 days). In collaboration with the Motus Wildlife Tracking System, we established an array of automated of sixteen radio telemetry towers positioned across a greater than150 km section of the U.S. Atlantic coast from Cape Cod, MA to Long Island, NY (see Fig. 0(a)). Each tower consisted of six Yagi antennas (Cushcraft PLC, 9element Yagi antenna), each having a horizontal plane beamwidth of 35 and whose main beams were separated in bearing by 60 and endmounted in a radial configuration atop a 12.2 m mast (see Fig. 0(b)). The antennas were connected via a coaxial cable to ports on a data logging receiving unit (SRX600, Lotek Wireless) that ran continuously using two 12V deepcycle batteries charged by a 120Watt solar system.
Each receiving unit was programmed to automatically log several types of data from each antenna, including: transmitter ID number, date, time (synchronized among all receivers in network using Global Positioning System (GPS) clocks), antenna (defined by receiving station and bearing), and power received (linear display scale: 0255). The received timeseries power signal was sampled sequentially by each Yagi antenna of a tower, each with a dwell time of 6.5 s (greater than the pulse duration). Thus the minimum separation between two successive readings at any Yagi antena is 32.5 s. However, the timeseries detection data can be nonuniformly spaced as tagged individuals traverse regions of the antenna beam from which signal strength is very low and not detected by the receivers.
We collected systematic calibration data of the reception pattern of the automated telemetry stations using a test transmitter attached to a partially frozen bird carcass that was suspended from a kite and flown from the back of a boat [25]. During our calibration surveys, the boat was traveled at a constant speed of approximately 10 m/s and the transmitter flew at an altitude of approximately 30 m.
To estimate the bird location based on these measurements, we develop a statespace based Kalman filtering technique, consisting of a system movement model and an observation model. The statespace modeling approach is similar to that employed by [26, 27], but the movement model was extended to accommodate correlated vertical motion based on a CoxIngersollRoss model [28] in addition to having correlated motion in the transverse plane, modeled by a modified OrnsteinUhlenbeck velocity model [29]. Conditions required for the dynamical system to be observable (uniquely solvable) are discussed. The observation model is further tailored to a Yagi antenna pattern conforming to the manufacturer’s catalog data. We conducted calibrations of signal strength measurements of the receiver to develop a linearized observation model that facilitates estimation by standard Kalman filtering. We then test the validity of the Kalman filter model on a simulated trajectory, and lastly apply the algorithm to predict the trajectory of a VHFtagged shorebird tracked over several hours by an automated radio telemetry array along the U.S Atlantic coast.
Our paper is organized as follows. The radio telemetry system is described in §2.1 and challenges specific to measurement system are described in §2.2. The movement model in the transverse plane and vertical direction are first developed in §2.3. The Yagi antenna pattern used in the observation model is constructed in §2.4.1, and making use of calibrated kite measurements, the receiver component governing signal gain is developed in §2.4.2. The linearized observation model for estimation by standard Kalman filtering is then developed in §2.4.3 and the Kalman filter algorithm is presented in §2.5. The validity of the models so developed is then tested in §3.1 on a controlled, synthetic bird trajectory. Finally, in §3.2, the algorithm is applied to predict the migratory bird trajectory using measurements gathered by our array of receivers.
2.2 Challenges and Observations
The sixteen towers that were erected in the Atlantic ocean region for the purpose of this study are shown in Fig. 0(a). Fig. 0(b) shows the typical set of 6 Yagi arrays erected at each tower. Given the relatively small size of the nanotag at the operating wavelength, it is reasonable to first assume that the radio nanotags have an omnidirectional pattern.
line in the left figure represents the axis of a Yagi array. The array axes are separated by
60 and the horizontal plane beamwidth is each array is .
Our measurement system and received data are characterized by the following challenges and observations which influence our model development:

Ch: Simultaneous measurements over multiple towers were not available, thus precluding triangulation.

Ch: The receiver records integer values in the range based on the power received on a decibel scale from proprietary coded signals.

Ch: The power received at any given time is not explicitly dependent on the instantaneous speed of nanotag. Thus bird speed cannot be directly inferred from power measurements alone.

Ob: When a signal is recorded by the receiver attached to a tower, it merely implies that the bird is within its range. Here, we are trying, in addition, to estimate the location of the bird in 3D space based on the intensity of the received signal. Assuming identical tags, the power received depends in a complicated way on the combined location coordinates of the moving transmitter and the radiation pattern of the Yagi array. The same signal power can be received from a multitude of locations even if the transmitter remains in the main beam of the receiving antenna. There is a good probability of even receiving the signals via the antenna side lobes and back lobe. Thus the inversion of three spatial coordinates from a single power measurement is necessarily nonunique.

Ob: Signal was received over highly nonuniform timeintervals. The minimum separation between successive measurements at any tower was around 5 s. The maximum separation between successive observations could involve multiple towers and exceeded over 12 hours in our case.

Ob: It was observed experimentally that during certain times the bird was detected via multiple beams and towers that were widely separated. This would be possible if the bird were flying at high altitudes. So a dynamic height dependence must be included in the theoretical models.
It might be tempting to consider a simple dynamic model consisting of three independent local level [30] systems along the three spatial coordinates. However, a quick analysis will reveal that it will render the system unobservable, thereby, resulting in highly illconditioned inversions. Observability means that the state vector can be uniquely determined from finite number of measurements [31]. One way to address the nonobservability quandary is to consider a firstorder system with decay and select varying decay constants in the threespatial directions. In addition to suffering from anisotropic behavior in space, such a first order system will result in very limited translational motion in space. We therefore consider a firstorder system of the OrnsteinUhlenbeck type in which it is the speed variable (and not the position variable) that driven by white noise. The corresponding displacement will be equal to the time integral of speed and will consist of a correlated random walk with a good component of drift. Our movement and observation models are characterized by the following specific features:

For the sake of analysis, the bird is treated like a particle that is described by its three spatial coordinates and any applicable speeds. The movement models in the transverse plane and in the vertical coordinate, , are different, somewhat mimicking the true flight pattern of a shorebird. In particular, the model permits nonrectilinear motion in the transverse plane and level motion with fluctuations in the vertical plane such that altitude always remains above the mean sea level. Additionally, model parameters governing flight can be optimized to reflect either directed or highly irregular movement.

In the transverse plane, the state space along each coordinate consists of location and speed as state variables. In the vertical coordinate the state space consists of location variable only.

The bird movement model has the required complexity to make all of state variables almost observable despite the receiver power being dependent only on the location coordinates as highlighted in observation . Making the system observable will also address the challenge by coupling the speed variables implicitly to power measurements.

The vertical coordinate comes into play both through the range (to tower) coordinate as well as through the multipath phenomenon. Presence of multipath will account for the height gain trend that is highlighted in observation . To make the dynamic model more physically reasonable we make the transverse coordinates and the vertical coordinate statistically correlated.

The overall system has the parameters tuned such that the model is moreorless isotropic in the transverse plane. The resulting equations in the final discrete time model are not subject to any approximations even for irregularly spaced measurements in time.

A receiver model is constructed that directly relates the received power to the display reading . This is based on nonlinear least square fitting of calibration data. This will address the challenge highlighted under . Saturation at the upper end caused by receiver nonlinearity and at the lower end by receiver noise is accounted for in the receiver model by the incorporation of a softlimiter and measurement noise respectively.

An analytical antenna pattern is synthesized that mimics the radiation pattern supplied by the manufacturer.

A nonlinear observation model is constructed from an analytical antenna pattern function, range, heightgain factor and the receiver model as determined in item (6) above.

Nonlinear Kalman filtering technique [32] is utilized to determine the mean trajectory of the bird based on the system model assumed and conditioned by the measurements.

When large time gaps in the received data are encountered, the model is initialized again with the initial coordinates for the new run set such that they are closest to the most recently determined mean coordinates. The overall mean trajectory is sensitive to the initial coordinates assigned to the model as well as to the strategy employed for handling large time gaps. In coming up with schemes for dealing with large time gaps, we keep in mind the absolute maximum speed that the bird can fly to narrow down on the many possibilities.
2.3 Correlated Motion Model
The dynamical variables in the transverse plane are the particle position, , and the corresponding instantaneous velocity components , whereas in the vertical plane the dynamical variable is such that the altitude coordinate is . We denote time derivative by an overdot on the variable. The equations of motion are governed by an OrnsteinUhlenbeck (OU) process [29], [33] without drift in the transverse plane and by the CoxIngersollRoss model [28] in the vertical direction:
(1)  
(2)  
(3)  
(4)  
(5) 
where are positive real constants with units of [], which determine the decay of , , and are white Gaussian noise processes[34] with the property that , being the Dirac delta function and likewise for , and , are positive constants with units [ms] that determine the variances of white noise and are real constants that determine the correlation between , and . It is assumed that the noise processes and driving the two Cartesian components are independent of each other and are also independent of . The operation denotes ensemble average and will involve all stochastic quantities. The first order equation (2) satisfied by the velocity variable with a decay constant and driven by random forces and is simply a statement of Newton’s force law. It is possible to incorporate environmental covariates such as wind speed and direction [26] into the present movement model by way of introducing potentials and the corresponding external forces^{1}^{1}1For a particle of mass with the potential energy, , the force per unit mass arising from the potential is . into the right hand side of (2). However, we will not consider external forces here. The first equation (1) relates the time rate of change of position to velocity. A fixed point of the system in the absence of noise is obtained by setting all time derivatives to zero. In the coordinate this results in for the model above, where is a constant. Fixed points are points of stagnation in the sense that if the state trajectory reaches it at some instant of time, it will remain there for all future times in the absence of noise. In an ordinary OU process a constant term is also added to the right hand side of the velocity equation (2). In that case the velocity settles around and the position will asymptotically reach in the absence of noise. Consequently, the transverse trajectories will tend be rectilinear in an ordinary OU process in the absence of noise. To remove this unwanted bias and to model nondirected movement, we assume here. For simplicity of notation, we will omit the argument from all quantities when there is no danger of confusion. Using Itô’s formula [35], [36] it can be shown that the altitude satisfies the nonlinear dynamical equation
(6)  
In mathematical finance, a model such as this is referred to as the CoxIngersollRoss model [28] that describes the evolution of interest rates (assuming they are ) over time. The fixed point for altitude is . Figure 2 shows a sample evolution of the altitude for ms, m and m. It is seen that the model permits finescale variation in altitude to capture the altitude variation of birds in flight. Because of the nonlinearity present in (6), fluctuations in the height variable will be more correlated than those of the intermediate variable . To maintain isotropy of coupling between to and to we choose .
Denoting the state vector , where a prime denotes matrix transpose and the noise vector by , the dynamical equations can be written in a matrix form as
(7) 
where the blockdiagonal system matrix is
(8) 
and the noise coefficient matrix is
(9) 
The parameters contained in the system are and are ten in number. Equation (7) represents a linear timeinvariant (LTI) system driven by random white Gaussian noise. Even though only the velocity components of the state vector and are directly driven by noise, the noise gets coupled to the other component by virtue of the system matrix being nondiagonal. However, because of linearity and timeinvariance of the system, it can be shown that the noise present in the components remains white Gaussian, albeit correlated. System (7) can be solved by employing techniques of stochastic calculus [36] or by using the traditional concept of transition matrix and Laplace transforms [37], [31]. Using any of these approaches, a recursive solution to (7) can be obtained as
(10) 
where
(11) 
where , , and is a zero mean, white Gaussian random vector with a covariance matrix :
(12) 
where is the Kronecker’s delta and is an identity matrix of size . Note that and as . It is also clear that is a Hermitian matrix, viz., , where superscript * denotes complex conjugation. The entries of the matrix can be easily found as the integrals involved are elementary. For instance,
(13) 
where . Equation (10) determines the state vector at time given the state vector and the stochastic input at time . It is an exact equation in the sense that it is not subject to any discretization error and does not require the time increment to be uniform. To gain some insight into the parameters and , let us consider the recursive solution of the velocity equation (2):
(14) 
where is a normal random variable.
For , and
(15) 
In this case, the updated velocity is relatively independent of the parameter and jumps from its previous value by a Gaussian random variable of standard deviation . This regime will lead to substantial displacements (integral of velocity) and may be a good choice for describing migration type of bahavior.
On the other hand for and , we have
(16) 
In this case the updated velocity is independent of the previous velocity and is mostly determined by random fluctuations and the parameter . As a result the displacement here will be small and this choice may be most suitable for describing nesting type of behavior. Here we may think of as the correlation time for the velocity component . For instance, when [s] the correlation time for velocity is about 1.1 hour.
It is desirable to keep the pairs and close to each other and choose appropriately so as to maintain isotropy of motion. To understand this, let the polar coordinates corresponding to be . Using , it is easy to see that
(18)  
which represents a nonlinear and a timevarying system. For any finitebandwidth physical system driven by white noise , the system response function is assumed to lag the white noise by an infinitesimal time interval [36]^{2}^{2}2This, indeed, forms the basis of Itô calculus.. Consequently, . Hence it is easy to see that , because is considered as a response function to and . For isotropy we require the variance of fluctuations in the radial velocity component to be independent of . It is easy to see that this will be accomplished when and . Fig. 3 shows a sample trajectory in the transverse plane with parameters as shown in the figure. Time is indicated as color on the trajectory. For the parameters chosen, the velocities remain in the range ms and show a behavior described by the regime in (15). It also seen that the position has a drift and is much more correlated in time than the velocity.
plane in the presence of system noise.
2.4 Observation Model
In this section we will first present a model that will express the power received in terms of the Cartesian coordinates of the tower and of the bird . It is known that the Lotek SRX600 receiving unit displays and outputs integer values that are proportional to the power received on a dBm (dBm = 10 power in milliwatts)) scale. Keeping this in mind we will then outline a procedure that will determine the constants involved and develop an observation model in terms of the state variables. Finally we will linearize the observation model that will permit application of the extended Kalman filter algorithm.
2.4.1 Yagi Array Pattern Model
For subsequent analysis, it is required to have knowledge of the radiation pattern of the Yagi antenna. Even though the radiation pattern values could be read off from the manufacturer data sheets, it is highly desirable to have an analytical form that will facilitate taking various gradients of the field strength later. The driven and the parasitic elements of the Yagi array all had a length, , of around half a wavelength (freespace wavelength m) and an overall array length of about 3.5 m. The array has an absolute gain of around 11.1 dB (decibels) and a fonttoback ratio of 20 dB [38]. If the main beam of the Yagi points in an azimuthal direction in a tower translated local coordinate system and are the polar coordinates of the point relative to the : , then the angle made by the point relative to the main beam is , Figure 3(a). The direct lineofsight distance between the tower and the bird is .
tal () plane.
For the purpose of developing the observation model, the array factor here is treated as arising from a continuous line source of effective length with uniform amplitude distribution and HansenWoodyard phase distribution [39]. In the plane, the component of the farzone radiated electric field is given by
(19)  
(20) 
where is a constant, , , , , and the phase per unit length for the HansenWoodyard condition. Choosing m yields a beamwidth of 35.8 and a frontback ratio of 22 dB, which are close to the values specified in the data sheets. The normalized field pattern versus is shown as solid line in Fig. 3(b) and closely matches the manufacturer’s pattern data despite the increased model length of m versus the actual array physical length of m. For subsequent analysis we write the power gain of the Yagi antenna as , where is some constant.
In the presence of multipath, equation (19) must be modified. If the signal arrives via two paths–one along a direct path between the transmitter and receiver having a path length , and an equally strong second one via ground reflection having a path length and reflection coefficient –then the received electric field is [40]
(21) 
It is assumed in the above equation that both rays arrive in a direction close to the main beam in the elevation plane so that all of the dependence is contained in and only.
2.4.2 Receiver Model
It was observed experimentally that during certain times the bird was detected via multiple beams and towers that were widely separated. This would be possible if the bird flies at high altitudes where there is substantial height gain associated with large ranges. A height gain will be apparent when there is (i) a lineofsight signal between the tower and the bird and (ii) a ground reflected signal arriving from low grazing angles [40] with . The appropriate expression for the received electric field in this situation is (21) which will result in a height gain factor. If the gain and transmitted power of the radio tag are and , respectively, then the received signal using the two ray model for horizontal polarizations in the absence of noise is [40]
(22)  
where that depends on the radio tag EIRP (Effective Isotropic Radiated Power ) and the receiving antenna gain factor . The factor in (22) involving the sine function is the height gain factor, which increases monotonically with bird height until the argument becomes equal to . For most cases encountered in practice the argument of the sine function will be less than . As an example if km, m, m the argument is 0.784. For a given range, the height at which the argument reaches is equal to . When km and receiver height m, m.
As stated previously, the SRX 600 receiver displays the signal strength, , using onebyte (i.e., displays integers ) with a scale monotonically proportional to the logarithm of the power received. At the low end, the received signal will be corrupted by receiver noise which will cause the display value not go down below certain minimum . At the upper end the receiver is designed to saturate causing the display not to exceed a maximum value, . The range of recorded display values will accordingly fall within the range . If the minimum receiver power (equal to the receiver noise power) is , the mean received power in the presence of noise is . We model the receiver nonlinearity by employing a soft limiter and express the limiter modified measurement, , as
(23)  
where the function is included for softlimiting, is the constant of proportionality and . Equation (23) may be rewritten in a simplified form as
(24) 
The inverse relation between and is easily obtained from (24) as
(25) 
Given reference measurements with known locations of the radio tag, the constants and may be estimated by performing nonlinear least square fitting. For a given the constant is equal to
(26) 
The constant is then found by minimizing the mean square error between the model prediction and actual data. We conducted reference measurements using a radio tag attached to a kite that was, in turn, attached to a boat whose coordinates were recorded using GPS signals. The boat was made to travel along a straight line in the direction of the main beam of a known tower antenna (MNYN, Fig.1a). The kite heights were approximately 100 ft and 200 ft along the two tracks. Signal from the 100 ft kite was measured by Southwest beam (beam5) and that from the 200 ft kite was measured by the Northwest beam (beam6). From the recorder time series signals we estimate the constants as . The value for agrees with the estimate based on thermal noise present in the receiver, where is the Boltzmann’s constant, is room temperature, is receiver bandwidth, and is its noise figure. For a radio tag with an EIRP of mW and a Yagi array with a gain of , we get when K, MHz, and . When , equation (24) gives , which for yields . So when the received power equals the noise power generated within the receiver, the display on Lotek SRX600 receiver is predicted to be 52.5. For received signal levels below the noise level, the display values below 52.5 will be unsteady and corrupted by noise.
2.4.3 Linearalized Observation Model
The instantaneous received power in the presence of noise is represented as
(27) 
where is a normal random variable that is statistically independent of the measurement . Clearly, the instantaneous power depends on the location coordinates and noise in a highly nonlinear fashion. The mean received power is equal to . The variance of the received power is . The mean and the variance are both seen to depend nonlinearly on the state vector and it is necessary to linearize the observation equation about a nominal trajectory before applying the standard Kalman filtering^{3}^{3}3The standard Kalman filter algorithm applies when both the system model and the observation model are linear in the statespace vector . However, when either one of the system or observation models is nonlinear in , the model must first be linearized about a nominal trajectory and then the linear algorithm can be applied. This is known as the extended Kalman filter. technique. Recall that
(28)  
(29)  
(30) 
Expanding in a Taylor’s series about a nominal trajectory and keeping only the linear terms in and noise (which is now labeled as ) yields
(31) 
where , , , , being the normal random variable, is the measurement row vector equal to
(32) 
2.5 Location Tracking by Kalman Filtering
Given the system model, (10), the observation model (31), a set of measurements , system and observation noise covariance parameters , , and , the conditional mean of bird location can be estimated by using the extended Kalman filtering technique as it evolves in time. What are recorded by the measurement system are the corrupted (by measurement noise) signal powers associated with a particular bird trajectory. The estimated mean bird trajectory will be dictated both by the system model and the measurements available. As a byproduct the variances associated with the mean trajectory are also calculated by the Kalman filter. For LTI systems operating in the presence of Gaussian noise, the Kalman filter gives the best estimate based on the current measurement. The following notation is first introduced:
covariance of measurement error  
The nominal trajectory for linearizing the measurement equation is chosen to be so that , . The discretetime extended Kalman filter, [32], [30] is given by the following equations, which are computed for each time step :
(33)  
(34)  
(35)  
(36)  
(37)  
(38)  
(39) 
It is important to note that when the measurement noise is very high, and the norm of the Kalman gain vector . Consequently, , meaning that the mean trajectory is governed primarily by the ‘predictor’ part of the Kalman filter or, equivalently, by the noisefree system model (33). On the other hand, when the system noise is very high, , and the ‘corrector’ part of the Kalman filter (38) will also gain prominence.
For the Kalman filter to estimate the states effectively, the system must be observable. Observability depends on the system coefficient matrix and the measurement matrix . For an LTI system to be completely observable, the observability matrix must be full rank [41]. However, given that the system models along the three spatial coordinates are uncoupled and power measurements are independent of the individual’s speed, the observability matrix can at best have a rank of 4 in our case. So that the system is almost observable (with a degree of observability of ) we have to choose the system parameters such that rank of does not fall below 4. A necessary condition for the rank to be at least 4 with is that . Hence isotropy must be sacrificed in order to meet the observability condition. We choose to meet the almost observability condition.
The more complete dynamic prediction that is achieved through the proposed dynamical model and Kalman filtering may be contrasted with the static prediction of the triplet by directly inverting equation (22) without recourse to any dynamic system. For a given value of measured signal , we have from (22) in a noisefree environment that
(40) 
It is apparent from (40) that a multitude of combinations can give rise to the same value of (and hence ) as we have observed in . For bird and tower heights much smaller than the horizontal bird range to tower, , we have the approximation and in (40), thus resulting in
(41) 
Hence the range squared determined for any given measurement follows the shape of the field radiation pattern of the Yagi array. The dotdashed line in Fig. 3(b) shows normalized to its maximum possible value, for a given measurement and bird height . Thus a given measurement value could arise from a nanotag located at large ranges within the 3 dB beam width of the array or from relatively shorter ranges located along the side lobes or back lobe of the array. Note also that the range scales as . So the same signal can be received from a bird that is, say, at 6 km along the main beam or at 1.1 km along a dB sidelobe or at 600 m along a dB back lobe. These numbers are representative of our Yagi array whose pattern is shown in Fig. 3(b). However, if the height and bearing of the bird are both known by some other means, then equation (41) can be used to uniquely determine the range to a measuring tower. Equation (41) also shows that for a given measurement , longer ranges to the tower are possible if the bird flies at higher altitudes as we have indicated in .
Under the assumption that the signals arrive primarily through the main beam of the receive array so that , we get
(42) 
where is the halfpower angle of the receiving antenna. (For the Yagi antenna used in our measurements .). Equation (42) constitutes what we label as the static model.
3 Results
3.1 Validation Using Simulated Bird Trajectory
Assuming tower coordinates of and (coordinate system UTM Zone 19N, units in m) we generated a simulated trajectory using equations (10) for a bird starting at with velocity components , and an initial height of m (equal to the tower height above sea level). Locations were simulated every 6 seconds (equivalent to the burst rate interval of the nanotag transmitter) for a total of 1,200 seconds. Then with the assumed tower coordinates and choosing the appropriate beam, we calculate the horizontal distance, height and bearing. The corresponding values of and are calculated using (22), (23) and (24). In the Kalman filter estimation, we have excluded points that produce the display values which correspond to received signaltonoise ratio of less than dB. The selected values are then input to the Kalman filter and state vector estimated for each instant of time. Fig 5 shows the comparison between the estimated trajectory and the actual trajectory in both planes. It is to be noted that the conditional mean trajectory estimated by the Kalman filter is not expected to exactly follow the actual trajectory, which is all but one realization of a whole ensemble of possible trajectories. However, the predicted trajectory is seen to capture the long term trends quite well including turnings in the horizontal plane even though the range at end of the interval falls short by about 1.6 km. For the predicted trajectory we calculate the display values as computed from (24) and compare them to the actual Zvalues in Fig.6. It is again seen that the display values follow the trends very well.
ory. Distances shown are relative to tower coordinates . The various paramet
ers used in the simulation were ,
.
3.2 Estimated Trajectory for a Migrating Piping Plover
We next apply the algorithm to predict the location of an adult Piping Plover (Charadrius melodus) that was tracked using multiple beams/towers during its migratory departure from the breeding grounds Fig. 0(a). The Piping Plover was captured during the incubation period using a walkin trap [42] at its nest on Aquidneck Island, Rhode Island, USA. The plover was fitted with a nanotag that was attached to its interscapular region using epoxy. Signals from the bird arrived at the towers at highly irregular intervals interspersed with several large gaps. We will show the predictions carried out over one large subinterval during which the bird was thought migrating in a southwest direction. The radio tag on the bird was detected by several towers that included SACH, TRUS, NAPA, BISE, PLIS, and MNTK (see Fig. 0(a)). Several hours later it was detected in New Jersey by a tracking station on the Motus network (www.motus.org). As remarked earlier, numerous trajectories can all lead to the same set of measurements. Because tracking is recursive in nature and is intimately tied to the initial state, some strategy must be employed to narrow down on the number of possible trajectories. Different strategies will lead to different estimated trajectories because of the nonuniqueness of the problem at hand. In the results we show, we have employed the following assumptions and strategies:

An initial height is chosen based on previous or expert knowledge.

The system parameters are chosen based on a trail and error basis. From isotropy and observability criteria . Our numerical trials have indicated that a good rule of thumb for modeling migratory type of behavior is . From (15) it is seen that the condition is satisfied as long as the maximum time gap between successive data arrivals at an antenna satisfies s. For a reasonable random change in speed of about 1 m/s between successive data observations in the transverse plane, the noise parameters may be chosen according to following (15). For example with , s, one gets the values [ms] and [ms]. Similar order of values for , and also implied from (17).

A maximum possible radial speed for the bird is specified.

An initial covariance matrix is chosen. A reasonable choice is to choose to be diagonal with entries proportional to the uncertainties of the state vector. For example, we chose here.

The initial state is determined by considering the first two measurement points recorded at times