Monte Carlo Simulation Analysis of Tagged Fish Radio ... - Springer Link
J Intell Robot Syst (2014) 74:287–307 DOI 10.1007/s10846-013-9949-9
Monte Carlo Simulation Analysis of Tagged Fish Radio Tracking Performance by Swarming Unmanned Aerial Vehicles in Fractional Order Potential Fields Austin M. Jensen · David K. Geller · YangQuan Chen
Received: 1 September 2013 / Accepted: 13 September 2013 / Published online: 18 October 2013 © Springer Science+Business Media Dordrecht 2013
Abstract Tracking fish using implanted radio transmitters is an important part of studying and preserving native fish species. However, conventional methods for locating the fish after they are tagged can be time consuming and costly. Unmanned Aerial Vehicles (UAV) have been used in general radio localization applications and can possibly be used to locate fish quickly and effectively. However, the methods developed for multi-UAV navigation and transmitter localization are complex and might not work well for practical and routine use. This work focuses on developing simple methods for multi-UAV navigation and transmitter localization. A real-
A. M. Jensen (B) Utah Water Research Laboratory, Utah State University, Logan, UT 84321, USA e-mail: [email protected] URL: http://aggieair.usu.edu D. K. Geller Mechanical and Aerospace Engineering Department, Utah State University, Logan, UT 84321, USA e-mail: [email protected] Y. Chen Mechatronics, Embedded Systems and Automation (MESA) Lab, University of California, Merced, Merced, CA 95343, USA e-mail: [email protected]
world simulator is created to test these methods; it includes a signal propagation model based on actual data from a UAV. Swarm-like navigation methods (using potential fields) are used for multi-UAV navigation, and an Extended Kalman Filter is used, along with a simplified version of the propagation model, to estimate the location of the transmitter. Multiple navigation methods are introduced and compared using Monte Carlo Analysis. Despite a noisy signal and a simplified measurement model, the different navigation methods are able to estimate the location of the transmitter with one or more UAVs. Keywords Unmanned aerial vehicle · UAV · Swarm · AggieAir · Potential field navigation · Kalman filter · Radio localization · Fish tracking · Monte Carlo analysis · UAV simulation · Wireless propagation model
1 Introduction In order to preserve our native fish species in the watersheds of the United States, managing, monitoring, and studying fish habitat is very important. Extensive fish habitat studies are conducted to minimize environmental impact before building a dam, while restoring natural habitat in a river, or in ensuring the economic strength of a fisheries market. An important part of these studies is
288
tracking where these fish are migrating, moving and living. This is commonly accomplished by tagging the fish of interest with a radio transmitter then using the transmitter to locate the fish in the future to see where it has moved [4, 5]. Tagged fish are located using a receiver and a directional antenna. As the antenna is pointed in the direction of the fish, the signal sent from the transmitter becomes louder. The operator can hear the beacon through the receiver and navigates toward the fish. While this method has proven to be effective, it can be time consuming and costly. An alternative method for tracking fish utilizes an unmanned aerial vehicle (UAV) developed at Utah State University (USU) called AggieAir [7]. Unlike most UAVs, AggieAir is primarily a remote sensing platform. It carries many types of payload, including a visual camera, a near-infrared (NIR) camera, a thermal infrared camera [8], and air quality sensors. A fish tracking payload has also been recently developed for AggieAir [11]. Using the fish tracking payload, AggieAir is able to fly over a body of water and receive the signal transmitted by the tagged fish. As AggieAir flies over the area of interest, the fish tracker measures the signal strength from its receiver and logs other data, including position, orientation, and speed of the UAV. With this information, the UAV can be directed toward areas of the map to improve the estimate of the tag location. Others have looked into similar scenarios. In a general radio tracking problem, Frew et al. [2] proposed a radio localization technique cast as a distributed estimation problem. Consensus within a group of UAVs was reached with a decentralized method; however, the algorithms used are computationally complex and only considered communication with one neighboring UAV. Korner et al. [10] introduced a method that uses particle filters to track multiple targets with radio tags. Implementation employed only one UAV. While these methods work, a simple, decentralized method is needed for multi-UAV navigation and radio tag localization. A method employing these features could be more practical for use in small UAVs with embedded systems. This paper is an extension of the author’s previous work [6]. This work is focused on simple methods for UAV navigation, multi-UAV
J Intell Robot Syst (2014) 74:287–307
coordination, and radio transmitter location estimation. Simplicity is important for real-world, real-time fish tracking applications. Therefore, swarm-like rules (using potential fields [3]) were established for navigation and multi-UAV coordination/avoidance. Much has been added in this paper. Real data was collected from a UAV and used to develop a propagation model. This propagation model was used to develop a real-world simulation to generate signals for the UAVs. A simplified version of the propagation model and an Extended Kalman Filter were used to estimate the position of the transmitter. Four simple navigation methods to control one or more UAVs are also introduced. These methods generate potential fields based on the estimated location of the transmitter and the covariance matrix of the estimated position. Comparisons are made between the methods by testing them with the simulator and Monte Carlo analysis. The results give clear answers to which methods are effective, and what the expected performance would be in real-world scenarios.
2 Methods In order to properly test different methods to locate tagged fish, a real-world simulation was created. Figure 1 shows a diagram that describes the different pieces of the simulation. Multiple UAVs can be included in the simulation. Each UAV (i) uses wind (|W|, ∠W) and noise (ω, η) to generate ˆ i ). The truth states (Xi ), and estimation states ( X potential field map receives the estimation states from each of the UAVs and updates the map with their current position. If the UAV receives a
Fig. 1 Flight simulation system diagram
J Intell Robot Syst (2014) 74:287–307
289
signal from the fish tag, the UAV estimation states are also used to update the estimated tag position. The truth state of the UAV is used by the tag propagation model to determine whether or not a UAV is close enough to receive a measurement and what that measurement will be ( S˜ i ). When the ˆ T ) is updated with a new estimated tag state ( X measurement, the estimated tag state and the corresponding covariance matrix (PT ) are both used to update the potential field map. The potential field map provides a force (F) which the UAVs will use to control their heading. The UAVs share each other’s position, along with the estimated position and covariance matrix of the tag, indirectly through the potential map. Static objects can also be contained in the potential map like ground based obstacles and a boundary. The tag state estimation error (δ XT ) is used to determine the success of the tag estimation and is calculated by subtracting the tag state estimate with the tag truth state (XT ). Monte Carlo analysis was used with δ XT to determine whether or not a method was successful. The basic idea of Monte Carlo analysis is to run a high-fidelity simulation many times. Because of the noise and disturbances included in the simulation, each trajectory generated the simulation will be different. Monte Carlo analysis is used to show how the trajectories vary and predict expected performance in real-world scenarios. For example, the simulation in Fig. 1 was run 100 times to look at how much δ XT varies for the different methods. One single run by this simulation will be referred to as a flight. Equation 1 shows how an error trajectory δ Xk (t) is calculated for flight number k:
The mean of the error trajectories describes the accuracy that most of the flights should achieve for that given method. The variance describes the best and the worst that can be expected around the mean. The uncertainty of the variance is also considered. Many applications require thousands of runs to reduce the uncertainty of the variance. For this application, 100 flights is enough to compare navigation methods. In addition, the uncertainty of the variance is documented and confidence intervals around the variance curves are calculated using a chi-squared distribution. The subsections below describe the important pieces of the flight simulation. These pieces include the environmental dynamics, the UAV system dynamics, simulating and estimating the tag location, and UAV navigation using potential fields.
ˆ Tk (t). δ Xk (t) = XT (t) − X
b ∠W b˙ ∠W = − + ω∠ W . τ∠W
(1)
After all of the flights are finished, the error trajectories are compiled in (t): (t) = [δ X1 (t), δ X2 (t) · · · δ X N (t)]T .
(2)
Then the mean μ(t) and variance (t)2 of (t) are found for each time (t): μ(t) = E[(t)],
(t)2 = E[(μ(t) − (t))(μ(t) − (t))T ].
2.1 Environmental Dynamics Equation 5 is used to calculate the magnitude and direction of wind. |W|0 and ∠W0 are nominal values of the wind and b |W| and b ∠W are random biases modeled as first order Markov stochastic processes (Eq. 6). |W| = |W|0 + b |W| , ∠ W = ∠ W0 + b ∠ W ,
(5)
˙ 0 = 0, |W| ˙ 0 = 0, ∠W b |W| b˙ |W| = − + ω|W| , τ|W|
The driving noise (ω|W| , ω∠W ) is Gaussian with the following properties. E[ω|W| ] = 0,
E[ω|W| ω|W| ]=
E[ω∠W ] = 0,
E[ω∠W ω∠ W] =
(3)
(4)
(6)
2 2σ|W|
τ|W|
δ(t − t ).
2σ∠2 W δ(t − t ). τ∠W (7)
290
J Intell Robot Syst (2014) 74:287–307
2.2 UAV System Dynamics Figure 2 shows a detailed diagram of the simulated UAV system for U AVi . The force from the potential field map (F) is used by the controller to generate ui (desired rate of change for the heading), which is used along with noise and wind to generate the UAV truth states (Xi ) and GPS measurements (z˜ i ). The GPS measurements are used as measurements in the Kalman Filter to generate the estimated states of the UAV. 2.2.1 Truth Dynamics
Equation 10 shows the equations to calculate other important truth variables: ground speed (vgi ), and heading (φi ). vgi = x˙ 2i + y˙ 2i , φi = tan−1
y˙ i . x˙ i
(10)
Part of the truth model includes generating the GPS measurements that will be used by the navigation. As shown by Eq. 11, the kth measurement is a function of the current position (xi or yi ) and the respective random bias: x˜ i [k] = xi + b xi ,
Equation 8 shows the UAV truth state dynamics. xi and yi describe the position of the UAV. The time derivative of the position is a function of the airspeed vi , the yaw ψi , and the x and y components of the wind (Wx , W y ). The airspeed is assumed to be constant and the change in yaw is equal to the controller output (ui ). x˙ i y˙ i v˙ i ψ˙ i b˙ xi b˙ yi
= = = =
vi cos(ψi ) + Wx , vi sin(ψi ) + W y , 0, ui , b xi =− + ωxi , τx b yi =− + ω yi . τy
E[ω yi ] = 0,
2σx2 δ(t − t ). τx 2σ y2 E[ω yi ωyi ] = δ(t − t ). τy
(11)
A measured ground speed (v˜ gi [k]) and heading (φ˜i [k]) are also generated from the GPS measurements: (x˜ i [k]− x˜ i [k−1])2 +( y˜ i [k]− y˜ i [k−1])2 v˜ gi [k] = , t[k]−t[k−1] φ˜ i [k] = tan−1
y˜ i [k] − y˜ i [k − 1] . x˜ i [k] − x˜ i [k − 1]
(12)
2.2.2 UAV Navigation Estimation
(8)
b xi and b yi are biases in the GPS measurement. These biases are also modeled as first order Markov stochastic processes. E[ωxi ] = 0,
y˜ i [k] = yi + b yi .
E[ωxi ωxi ] =
To control the heading of the UAV according to the potential field map, an estimated position ˆ This is and heading of the UAV is needed (X). achieved through a Kalman filter. ˆ = [xˆ i yˆ i φˆi ]T . X
(13)
To develop the Kalman filter, a design model is used where Xdm = [xdm ydm φidm ]T , i i
(14)
and (9)
dm x˙ dm = v˜ gi cos(φidm ) + ωdm i x , dm y˙ dm = v˜ gi sin(φidm ) + ωdm i y ,
φ˙ idm = ui + ωφdm .
(15)
dm dm ωdm x , ω y , and ωφ represent the process noise of the propagation. Equation 15 can also be represented in the following vector form:
Fig. 2 UAV system diagram
dm ˙ dm = fˆ(Xdm , v˜ gi ˆ dm , X , ui ) + Bw
(16)
J Intell Robot Syst (2014) 74:287–307
where fˆ =
291
where
⎡
⎤ dm v˜ gi cos(φidm ) dm ⎣ v˜ gi sin(φidm ) ⎦ ,
ˆ X ˆ − ), zˆ i [k] = h( ⎡ ⎤ 100 ˆ ∂h Hˆ = = ⎣0 1 0 ⎦ , ˆ ∂X 001
ui ⎤
⎡
100 Bˆ = ⎣0 1 0⎦ , 001 dm
w
=
dm dm T [ωdm x ω y ωφ ] .
(17)
(26) (27) ⎡
⎤ σx2 0 0 Rˆ = E[ν dm (ν dm )T ] = ⎣ 0 σ y2 0 ⎦ . 0 0 σφ2
(28)
Using this design, the filter state is propagated according to the following equations:
2.2.3 Heading Control
˙ˆ ˆ v˜ gi , ui ), X = fˆ(X,
(18)
˙ ˆ Bˆ T , Pˆ = Fˆ Pˆ + Pˆ Fˆ T + Bˆ Q
(19)
All of the UAVs share their location to the potential field map. When U AVi shares its location (xi and yi ), the potential field map finds the gradient of the potential in that area (created from other UAVs, the border, other obstacles, etc) and returns a force Fi (xi , yi ). The x and y components of the force are used to find the desired heading φ f i of the UAV (Eq. 29).
where
⎡ ⎤ 0 0 − sin(φˆ i ) ˆ ∂f Fˆ = = ⎣0 0 cos(φˆ i ) ⎦ , ˆ ∂X 00 0 ⎡ 2 ⎤ σ px 0 0 2 ˆ = E[wdm (wdm )T ] = ⎣ 0 σ py 0 ⎦. Q 2 0 0 σ pψ
(20)
(21)
For the update step of the Kalman filter, Eq. 22 shows how the GPS measurements are related to the states of the design model. dm x˜ dm + ηbdmx , i [k] = xi dm y˜ dm + ηbdmy , i [k] = yi
φ˜ idm [k] = φidm + ηφdm .
(22)
The vector form of Eq. 22 is shown below: ˆ dm ) + ν dm , z˜ i [k] = h(X where
⎡
⎤ dm
xi ⎥ ˆ dm ) = ⎢ h(X ⎣ ydm i ⎦, φidm
(23) ⎡
⎤ dm
ηb x ⎢ ⎥ ν dm = ⎣ηbdmy ⎦ . ηφdm
(24)
Fi (xi , yi ) = Fx ˆi + F y ˆj,
Fy φ f i = tan . Fx
Figure 3 shows a diagram of the closed-loop heading controller. The error (φe = φ f i − φˆ i ) is multiplied by a control gain k in order to give the appropriate amount of control to turn the UAV in the right direction. This control output (ui ) is the turning rate of the UAV and is limited to physical constraints (|ui | ≤ umax ). This physical constraint is applied to the control output before it is given to the UAV. 2.3 Simulating and Estimating the Tag Location ˜ received by the UAVs To simulate the signal ( S) from the fish tag, actual flight data was collected with a UAV to find the propagation model. This data was fit to Eq. 30 where fs is the signal model, r is the distance from the UAV to the tag (xt ,yt ),
Based on this design, the update equations for the Kalman filter are shown below: ˆ −1 , K = P− Hˆ T ( Hˆ P− Hˆ T + R) ˆ − (I − K H) ˆ T + K RK ˆ T, P+ = (I − K H)P ˆ+ =X ˆ − + K(z˜ − zˆ ), X
(25)
(29)
Fig. 3 Closed loop heading controller
292
J Intell Robot Syst (2014) 74:287–307
α is the difference between the polarization angles of the tag and the UAV, and ηs is Gaussian white noise. As Eq. 33 shows, α is symmetric around its own axis and the axis perpendicular to it (0 ≤ α ≤ π/2). S˜ = fs (r, α) + ηs .
(30)
E[ηs ] = 0,
(31)
E[ηsi ηsjT ] = σs2 δ(ij).
r = sqrt(xˆ t − xˆ )2 + ( yˆ t − yˆ )2 ),
(32)
⎧ ⎪ αu − αt + 180 ⎪ ⎪ ⎪ ⎨α − α t u α= ⎪ α − α u t ⎪ ⎪ ⎪ ⎩α − α + 180 t u
.
2.4 Navigation Using Potential Fields 2.4.1 Conventional Potential Fields Potential fields (PFs) are commonly used to dynamically plan the motion of ground [3] and air-based [13] mobile robots. They are used to attract or repel agents to or from objects and other agents. For example, if each agent was surrounded by a repulsive PFs they would all repel each other and avoid collision. Equation 39 shows a conventional definition of an attractive PF [9].
−180 ≤ αt − αu < −90 −90 ≤ αt − αu < 0 0 ≤ αt − αu < 90 90 ≤ αt − αu ≤ 180
(33) An Extended Kalman Filter was used to estimate the position of the tag. Xt is the state vector for the filter. x Xt = t . (34) yt ˙ t = 0); It is assumed that the tag will not move (X therefore, there is no need for a propagation step with this filter. However, a linearized measurement model is needed. Equation 35 shows the measurement model that was used in the Kalman Filter. hˆ is a simplified version of fs and is not a function of α. ˆ Sˆ = h(r).
(35)
With H, the Kalman gain, K, can be computed, and the state and covariance matrices can be updated using Eq. 36: K = P− H T (H P− H T + R)−1 , P+ = (I − KH)P− (I − KH)T + K RK T , ˆ ˆ+ =X ˆ − + K( S˜ − S), X t t
(36)
where, ∂h H= , ∂Xt
(37)
R = E[ηs (ηs )T ] = σs2 .
(38)
U(x) =
1 k(xa − xg )2 . 2
(39)
This PF is applied to the agents through a force defined as the following negative gradient of the PF: F(x) = −∇U(x).
(40)
As xa (position of agent) gets further away from xg (position of a goal), U(x) goes up. The slope of U(x) also increases as U(x) goes up, which increases the force acting on the agent from the PF and attracts the agent toward xg . A repulsive PF would increase as the distance from the object to the agent decreases and would repel the agent. 2.4.2 Fractional Order Potential Fields The method used to generate the PFs for this project is called the Fractional Order Potential Field (FOPF) [1]. The main advantage of the FOPF over other methods is that the FOPF allows the user to design repulsive PFs based on a level of danger. For example, a collision with another UAV would be very dangerous and therefore would get a different PF that produced a higher force than the PF of a less dangerous object like a boundary. A repulsive FOPF is derived starting with the definition of the Coulombian electric field E(r) E(r) =
q . 4π0 r 2
(41)
A single integration of E(r) produces the Coulombian potential from the electric field of a
J Intell Robot Syst (2014) 74:287–307
293
punctual charge (V1 (r)), and a double integration produces the Coulombian PF generated by the uniformly distributed charge along a straight-line segment V2 (r): V1 (r) =
q , 4π0 r
(42)
q ln r V2 (r) = . 4π0
(43)
Both V1 (r) and V2 (r) could be used as repulsive PFs. Figure 4 (for n = 1 and n = 2) shows the normalized potential of V1 (r) and V2 (r) vs. distance and shows that V2 (r) will have a stronger force at a larger distance than V1 (r). For a more dangerous object, V2 (r) would be preferred over V1 (r). It turns out that with each successive integration of E(r), more force is applied at a further distance. Fractional calculus can be used to generate these integrals, as well as fractional order integrals. This conversion using the Weyl fractional integral [12] is given below: Vn (r) = W r E(r) =
q 4π0 (n)
r
∞
(θ − r)n−1 dθ, r2 (44)
Fig. 4 Repulsive Fractional Order Potential Field (1 ≤ n ≤ 5)
where Vn (r) is the nth integral of E(r) (n can be any number greater than 0) and (n) is the Gamma function. 1 n−1 1 (n) = ln dt. (45) t 0 After manipulating Eq. 44, Vn (r) can be written as shown below: Vn (r) =
q (2 − n) ∀n ∈ (0, 2) (2, ∞). 4π0 r 2−n
(46)
To further simplify, Vn (r) can be normalized between 0 and 1 with a maximum and minimum distance (rmax and rmin , respectively) where the normalized potential function U n (r) would be 1 at the minimum distance and 0 at the maximum. Urep (r) =
n−2 Vn (r) − Vn (rmax ) r n−2 − rmax = n−2 . n−2 Vn (rmin ) − Vn (rmax ) rmin − rmax
(47) Figure 4 shows the result of Eq. 47 with 1 ≤ n ≤ 5. When n = 2, the denominator of Eq. 47 goes to zero and produces a singularity; therefore, Eq. 43 is normalized and used in its place. As Fig. 4 shows, an FOPF with a greater n will be more
294
J Intell Robot Syst (2014) 74:287–307
useful in dangerous cases and will apply a greater force sooner than an FOPF with a smaller n. One other advantage of the FOPF over other conventional PFs is that they are normalized and have maximum range. This will keep all PFs equal and unable to affect the agent if it is outside its range. For these reasons FOPFs are also used as attractive PFs (Fig. 5). To find the attractive FOPFs, Eq. 47 is mirrored around the x-axis. U att (r) = −
n−2 r n−2 − rmax n−2 n−2 rmin − rmax
.
(48)
2.4.3 Navigation for Tag Estimation Besides avoiding other UAVs and staying within bounds, the main navigation objective is to fly in areas to improve the estimate of the tag position. Four methods are presented: The Attractive Method, the Repulsive Method, the Offset Repulsive Method, and the Dual-offset Repulsive Method. Figure 6 shows a PF map for the Attractive Navigation method. In this method, a large attractive PF is formed at the estimated position of
Fig. 5 Attractive Fractional Order Potential Field (1 ≤ n ≤ 5)
the tag. This draws the UAVs to the estimated position (the center of the graph in Fig. 6). The Repulsive Method also has a large attractive field centered at the estimated position of the tag. However, it includes a small repulsive PF, also centered at the estimated position of the tag. The attractive field draws the UAVs in, and the repulsive field keeps them from getting close to the estimated position. Instead of going to the estimated position, they circle around it. Figure 7 shows the PF map for the Repulsive Method. The wide band 400 m from the center of the map is where the repulsive field meets the attractive field. This is the area of least potential and the location to which the UAVs will be attracted. The Offset Repulsive Method is identical to the Repulsive Method, except the small repulsive field is centered around a point offset from the estimated position of the tag. Figure 8 shows an example of the Offset Repulsive Method. In this example, the repulsive PF is placed at an offset of 100 m in the positive x direction. This offset creates the moon shaped area displayed in Fig. 8 to which the UAVs are attracted. The axis on which the offset is placed is dependent on the covariance of the estimate. If the offset is currently placed
J Intell Robot Syst (2014) 74:287–307
295
Fig. 6 Potential field map for the Attractive Method
Single Attractive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
200
400
600
800
1000
x(m)
on the x-axis and the UAVs are flying in that area, then the variance of the estimate along the x-axis will decrease. Once it is 10% lower than the variance along the y-axis, the offset will move clockwise to the y-axis and the UAVs will move with it. After a few iterations, the offset and the
Fig. 7 Potential field map for the Repulsive Method
UAVs will have made a full circle around the estimate in a clock-wise direction. The Dual-offset Repulsive Method uses two repulsive PFs. Each one is offset an equal distance from the estimated position of the tag to create the map shown in Fig. 9. In this specific
Single Attractive with Single Repulsive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
x(m)
200
400
600
800
1000
296
J Intell Robot Syst (2014) 74:287–307
Fig. 8 Potential field map for the Offset Repulsive Method
Single Attractive with Offset Repulsive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
200
400
600
800
1000
x(m)
example, the two offsets are placed 200 m away from the estimate on the x-axis. This creates the two low-potential areas shown above and below the repulsive PFs on the y-axis. Like the Offset Repulsive Method, this method also chooses the axis to place the offsets based on the variance of the estimation. If the variance along the x-axis is too high, then the offsets will be placed on the y-axis to attract the UAVs to the x-axis. This method might produce better results with multiple UAVs than the Offset Repulsive Method since
Fig. 9 Potential field map for the Dual-offset Repulsive Method
it encourages the UAVs to spread out and fly in different areas.
3 Results 3.1 Tag Signal Propagation Model To model the signal propagation, the AggieAir fish tracking payload was flown in the flower
Single Attractive with Dual Offset Repulsive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
x(m)
200
400
600
800
1000
J Intell Robot Syst (2014) 74:287–307
297
the signal strength (approx. 80 dB). This is to be expected in a wireless system, especially while using a relatively simple model. A more complex model may yield lower noise but is not needed for the simulation since the goal is to recreate the signal (even if it is very noisy) and estimate its location with a model that is even more simple. Once the measurement model for the signal was known, hˆ and H from Eqs. 35 and 37 could be specified. As shown in Eq. 50, only the first two terms in r were used for estimation; none of the α terms were used.
Fig. 10 The flower flight path used to collect flight data (footnote 1)
pattern shown in Fig. 101 at 400m above the tag. Before the flight, the position of the tag (xt , yt ) and the polarization angle (αt ) were both noted. During the flight, the payload software logged the signal strength and the position and orientation of the UAV. These data ware used to find r and α according to Eqs. 32 and 33. r and α were then used with the measured signal strength ( S˜ m ) to find the coefficients of fs (Eq. 49). fs (r, α) = β1 + β2r + β3r 2 + β4 r 3 + β5 + β6 α + β7 α 2 + β8 α 3 .
(49)
Figure 11 shows the result of the modeling. The surface in the plot is the model and the points around the surface are the measured signals received during the flight. Figure 12 shows the residuals between the model and the data. The data from this figure show the noise of the signal and are used to generate ηs in Eq. 30 for the simulation. Notice that the 3-sigma value for the noise is half of the maximum possible value for
1 Figure 10 does not show the path of the UAV during the simulations.
ˆ h(r) = β1 + β2r,
(50)
∂ hˆ β2 (xˆ t − xˆ ) β2 ( yˆ t − yˆ ) H= = , . ∂Xt r r
(51)
3.2 Flight Simulation Figure 13 shows the a map of the flight simulation environment. A border surrounds the flight area and has a potential field to keep the UAVs inside. The tag is located in the middle of the map (0,0) and the boundary of the signal is indicated by the black dotted circle around the tag location. The location of the UAVs are indicated on the map by a star symbol. The contours surrounding the UAV are their potential fields, and the blue line behind the UAV is its track. The estimated location of the tag is drawn on the map by the ‘o’ symbol. The large contour lines centered at the estimated tag location are the potential fields of the estimate (Single Attractive Method). The covariance of the tag location estimate is also indicated on the map by the covariance ellipse drawn in black. Initial simulations were run to validate the realworld disturbances added to the simulation (e.g. wind, GPS error, etc). Table 1 contains the values for the simulation parameters that were used for the following results. 3.2.1 Simulated Disturbance and Noise Results Figure 14 shows an example of the position error of the UAV in the x- and y-axis. The simulation did well at creating conditions that a UAV might encounter during flight.
298
J Intell Robot Syst (2014) 74:287–307
Fig. 11 The two-dimensional, fourth-order polynomial model plotted with the test flight data
Actual Signal Data and Model 120 100 Model
Signal
80
Actual Data
60 40 20
0 1000 800 600 r
400 200 0
The wind was also well simulated (Fig. 15). This is even more apparent when looking at the plots of ground speed (Fig. 16) and yaw vs. heading (Fig. 17). As is expected under windy conditions, the ground speed is highly variable even though the airspeed was constant. The difference between the values for yaw and heading (crab
Fig. 12 Residuals of the data vs the signal model
80
100
60
40
0
20
α
angle) also illustrates the success in modeling the wind. Figure 18 shows a plot of the measured signal ˜ vs. the signal strength without noise strength ( S) ( fs ). A noisy, variable signal is expected given the noisy data acquired and the dependence on polarization angle.
Residual error plot. 3−sigma: 40.3372 70 30 60 25
50
20
α
40
30
15
20
10
10
5
0
0
100
200
300
400 r
500
600
700
J Intell Robot Syst (2014) 74:287–307 Fig. 13 Flight simulation map
299 Border Potential Field
1000 800 600
Covariance Ellipse
UAV Location & Potential Field
400 200 0 −200 −400
Estimated Tag Location
Tag Location
−600 −800 −1000
−1000
3.2.2 PF Navigation Method Results The Attractive PF Method is demonstrated in Fig. 13. As expected, the UAV was drawn toward the estimated position of the tag due to the large attractive field centered at the estimated position. Figure 19 shows a flight using the Repulsive Method. With this method, the UAV was encourTable 1 Simulation parameters Symbol
Value
Description
dt GPSdt Sdt 3σx 3σ y 3σφ
0.1 1 4 5 5 2
τx τy |W|0 3σ|W| 3σ∠ W τ|W| τ∠ W umax k xT , yT αT 3σs
10 10 5 5 5 50 50 10 10 (0,0) 45 40.3
Simulation step size (s) GPS Measurement interval (s) Tag interval (s) 3-sigma value for GPS bias (m) 3-sigma value for GPS bias (m) 3-sigma value for heading process noise (degrees) Time constant for GPS bias in x (s) Time constant for GPS bias in y (s) Nominal wind magnitude (m/s) 3-sigma wind magnitude bias (m/s) 3-sigma wind direction bias (degrees) Time constant for wind magnitude (s) Time constant for wind direction (s) Max turning rate for UAVs (deg/s) Gain for heading controller Position of tag (m) Polarization angle of tag (degrees) 3-sigma value for signal strength noise (dB)
−500
0
500
1000
aged to circle around the estimated position; however during a few simulations, it would get turned around by the wind and blocked from covering the entire circle. The Offset Repulsive Method (Fig. 20) was more successful at directing the UAV completely around the estimate. In addition, it yields an interesting variance plot (Fig. 21) that shows “steps” of improvement in each axis when that particular axis is being worked on by the UAV. The Dual-Offset Repulsive Method yielded similar results to the Offset Repulsive Method. Figure 22 shows a flight with two UAVs. The Dual-Offset Repulsive Method was successful at separating the UAVs and directing them toward the two separate areas of low potential. 3.3 Monte Carlo Simulations Table 2 shows the results of the Monte Carlo Simulations. The mean, the 3-sigma values, and the 99 % chi-squared confidence interval for the 3-sigma values are listed for all four methods. Each method was also simulated with one, two and three UAVs. With the highest mean and 3-sigma values, the Attractive method had the worst performance and is not worth looking into further. Focusing on the other three methods and the simulations
300
J Intell Robot Syst (2014) 74:287–307
Fig. 14 Plot of UAV position error in X and Y
Position Error, var x: 2.6454 var y: 1.3555 6 5 4 3
Error (m)
2 1 0 −1 −2 −3 −4
10
20
30
40
50 60 Time(s)
70
80
90
100
70
80
90
100
70
80
90
100
Wind Magnitude (m/s)
Wind Magnitude 10 8 6 4
0
10
20
30
40
50 60 Time(s)
Wind Direction Wind Direction (deg)
Fig. 15 Plot of wind magnitude and direction
0
195 190 185 180 175 170
0
10
20
30
40
50 60 Time(s)
J Intell Robot Syst (2014) 74:287–307
301
Fig. 16 Plot of actual and measured ground speed
Ground Speed 22 Truth Measured
20
Speed(m/s)
18
16
14
12
10
8
Fig. 17 Plot of yaw, heading and wind direction
0
10
20
30
40
50 60 Time(s)
70
80
90
100
Yaw, Heading, Wind Direction Plot 180 160 140
yaw(psi) heading(phi) est heading(phi) Desired phi wind
Degrees
120 100 80 60 40 20 0 −20 5
10
15
20 25 Time(s)
30
35
40
302
J Intell Robot Syst (2014) 74:287–307
Fig. 18 Plot of measured and modeled signal strength
Measured vs Modeled Signal 90 80 70
Signal
60 50 40 30 20 10 0
0
50
with just one UAV, the following observations are noted. 1. The two offset methods have very similar results and have a lower mean value than the Repulsive Method.
Fig. 19 Repulsive method simulation results
100
150 Time(s)
200
250
300
2. The Repulsive Method may have a higher mean than the offset methods, but the upper 3-sigma value is lower. 3. The Repulsive Method is the only one that has a lower 3-sigma value greater than zero.
1000 800 600
wind
400 200 0 −200 −400 −600 −800 −1000
−1000
−500
0
500
1000
J Intell Robot Syst (2014) 74:287–307 Fig. 20 Offset repulsive method simulation results
303 1000 800 600
wind
400 200 0 −200 −400 −600 −800 −1000
−1000
−500
4. The Repulsive has the smallest gap between the mean and the 3-sigma values and the smallest Chi-2 spread.
The conclusion from these observations is that, while the offset methods may achieve higher
Fig. 21 Offset repulsive method tag position variance plot
0
500
1000
accuracy estimates, they are less precise and have higher uncertainty. The Repulsive Method may not be able to achieve errors as low as the Offset Methods; however, is much more reliable. The error trajectory plots of the Repulsive Method and the Offset Method also confirm this conclusion (Figs. 23 and 24).
Tag XY Estimation Variance 900 x y 800
Variance
700
600
500
400
300
200
0
50
100
150
200
250 300 Time(s)
350
400
450
500
304
J Intell Robot Syst (2014) 74:287–307
Fig. 22 Dual-offset repulsive method simulation results (2 UAVs)
1000 800 600
wind
400 200 0 −200 −400 −600 −800 −1000
−1000
The table also shows that, with the Attractive and Repulsive methods, having multiple UAVs improved the estimation performance. The opposite is true with the Offset and Dual-offset methods. For these methods, their accuracy and confidence were reduced as more UAVs were
Table 2 Monte Carlo simulation results for tag error dispersion and 100 flights
All of the results are in meters
−500
0
500
1000
added to the simulation. Even so, the general improvement shown in the table from multiple UAVs is not great. The greatest contribution that multiple UAVs made to the estimate is illustrated in the error trajectory plots. Figure 25 shows the trajectory plot for the Repulsive Method and
Method
# of UAVs
Final mean
Upper 3-sig
Lower 3-sig
Chi-2 spread
Att
1
274
671
0
149
Rep
1
251
343
159
35
Offset
1
160
408
0
93
Dual
1
160
360
0
76
Att
2
249
577
0
123
Rep
2
256
325
187
27
Offset
2
168
411
0
98
Dual
2
186
436
0
101
Att
3
234
572
0
127
Rep
3
255
318
192
27
Offset
3
186
443
0
111
Dual
3
190
471
0
108
J Intell Robot Syst (2014) 74:287–307
305
Fig. 23 MC simulation with 1 UAV using repulsive method
Tag Estimation Error Trajectories 1000 Number of Runs/Errors: 100/0 Final Mean: 251 3−Sig Values, Upper:343 Lower:159 Chi−2 Spread: 34.6246
Estimation Error (m)
800
Mean 3 Sigma 99% Chi2 Bounds
600
400
200
0
−200
−400
0
500
1000
1500
2000
2500
3000
Time (s)
Fig. 24 MC simulation with 1 UAV using offset repulsive method
Tag Estimation Error Trajectories 1000 Number of Runs/Errors: 100/0 Final Mean: 160 3−Sig Values, Upper:408 Lower:−88 Chi−2 Spread: 93.3655
Estimation Error (m)
800
Mean 3 Sigma 99% Chi2 Bounds
600
400
200
0
−200
−400
0
500
1000
1500
Time (s)
2000
2500
3000
306
J Intell Robot Syst (2014) 74:287–307
Fig. 25 MC simulation with 3 UAV using repulsive method
Tag Estimation Error Trajectories 1000 Number of Runs/Errors: 100/20 Final Mean: 255 3−Sig Values, Upper:318 Lower:192 Chi−2 Spread: 26.8739
Estimation Error (m)
800
Mean 3 Sigma 99% Chi2 Bounds
600
400
200
0
−200
−400
0
500
1000
1500
2000
2500
3000
Time (s)
three UAVs. Compared to the Repulsive plot with one UAV (Fig. 23), the plot with three UAVs converges ten times faster than the plot with one UAV. Even though using multiple UAVs did not improve the error in the estimate, they can be used to find the tag in a shorter time.
4 Conclusions A real-world simulation was created to compare the proposed navigation methods in the task of using UAVs to estimate the position of a small transmitter. Comparisons were also made between using single and multiple UAVs. The simulation proved successful at implementing realworld disturbances like GPS bias, a variable ground speed, and a crab angle. The estimation of the tag was also successful despite the noisy signal and the simplified measurement model. The Monte Carlo analysis showed that the Offset and Dual-offset Methods could estimate the tag with the highest accuracy; however the confidence on those estimates was very low. The Repulsive Method was not able to estimate with the same accuracy but was very consistent. The use of multiple UAVs did not improve the estimation accuracy; in some cases, it made the estimates less
accurate. However, multiple UAVs were able to find the tag ten times faster than only one. Future work includes a sensitivity analysis to figure out which error sources are major contributors to the tracking error. With this information, meaningful improvements can be made on the platform to increase the accuracy. Acknowledgements The authors would like to acknowledge Dr. Mac McKee of the Utah Water Research Laboratory. Thank you for your continued support. Gratitude also goes to Miguel Leonardo and Ivan Jimenez for their efforts and skills in developing the fish tracking payload. This work is supported by Utah Water Research Laboratory MLF 2006-2013.
References 1. Chen, Y., Wang, Z., Moore, K.: Optimal spraying control of a diffusion process using mobile actuator networks with fractional potential field based dynamic obstacle avoidance. In: Proc. of the IEEE International Conference on Networking, Sensing and Control, pp. 107–112. IEEE (2006) 2. Frew, E.W., Dixon, C., Argrow, B., Brown, T.: Radio source localization by a cooperating UAV team. In: Infotech@ Aerospace (2005) 3. Ge, S., Cui, Y.: Dynamic motion planning for mobile robots using potential field method. Auton. Robot. 13(3), 207–222 (2002)
J Intell Robot Syst (2014) 74:287–307 4. Gray, R.H., Haynes, J.M.: Spawning migration of adult chinook salmon (oncorhynchus tshawytscha) carrying external and internal radio transmitters. J. Fish. Res. Board Can. 36(9), 1060–1064 (1979) 5. Hillyard, R.W., Keeley, E.R.: Distribution and spawning migrations of fluvial bonneville cutthroat trout in the Bear River, Idaho. Tech. Rep., Idaho State University (2009) 6. Jensen, A., Chen, Y.: Tracking tagged fish with swarming unmanned aerial vehicles using fractional order potential fields and Kalman filtering. In: Proc. of the International Conference on Unmanned Aircraft Systems (ICUAS) (2013) 7. Jensen, A.M., Hardy, T., McKee, M., Chen, Y.Q.: Using a multispectral autonomous unmanned aerial remote sensing platform (AggieAir) for Riparian and Wetland applications. In: Proc. IEEE Int. Geoscience and Remote Sensing Symp. (IGARSS) (2011) 8. Jensen, A.M., McKee, M., Chen, Y.: Calibrating thermal imagery from an unmanned aerial system— AggieAir. In: Proc. IEEE Int. Geoscience and Remote Sensing Symp. (IGARSS) (2013)
307 9. Khatib, O.: Real-time obstacle avoidance for manipulators and mobile robots. In: Proc. of the IEEE International Conference on Robotics and Automation, vol. 2, pp. 500–505. Institute of Electrical and Electronics Engineers (1985) 10. Korner, F., Speck, R., Goktogan, A.H., Sukkarieh, S.: Autonomous airborne wildlife tracking using radio signal strength. In: Proc. of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 107–112. IEEE (2010) 11. Leonardo, M., Jensen, A., Coopmans, C., McKee, M., Chen, Y.: A miniature wildlife tracking UAV payload system using acoustic biotelemetry. In: Proc. of the ASME International Design Engineering Technical Conferences & Computers and Information in Engineering Conference (2013) 12. Miller, K.: The Weyl fractional calculus. In: Ross, B. (ed.) Fractional Calculus and its Applications. Lecture Notes in Mathematics, vol. 457, pp. 80–89. Springer Berlin Heidelberg (1975) 13. Paul, T., Krogstad, T.R., Gravdahl, J.T.: Modelling of UAV formation flight using 3D potential field. Simul. Model. Pract. Theory 16(9), 1453–1462 (2008)
Monte Carlo Simulation Analysis of Tagged Fish Radio Tracking Performance by Swarming Unmanned Aerial Vehicles in Fractional Order Potential Fields Austin M. Jensen · David K. Geller · YangQuan Chen
Received: 1 September 2013 / Accepted: 13 September 2013 / Published online: 18 October 2013 © Springer Science+Business Media Dordrecht 2013
Abstract Tracking fish using implanted radio transmitters is an important part of studying and preserving native fish species. However, conventional methods for locating the fish after they are tagged can be time consuming and costly. Unmanned Aerial Vehicles (UAV) have been used in general radio localization applications and can possibly be used to locate fish quickly and effectively. However, the methods developed for multi-UAV navigation and transmitter localization are complex and might not work well for practical and routine use. This work focuses on developing simple methods for multi-UAV navigation and transmitter localization. A real-
A. M. Jensen (B) Utah Water Research Laboratory, Utah State University, Logan, UT 84321, USA e-mail: [email protected] URL: http://aggieair.usu.edu D. K. Geller Mechanical and Aerospace Engineering Department, Utah State University, Logan, UT 84321, USA e-mail: [email protected] Y. Chen Mechatronics, Embedded Systems and Automation (MESA) Lab, University of California, Merced, Merced, CA 95343, USA e-mail: [email protected]
world simulator is created to test these methods; it includes a signal propagation model based on actual data from a UAV. Swarm-like navigation methods (using potential fields) are used for multi-UAV navigation, and an Extended Kalman Filter is used, along with a simplified version of the propagation model, to estimate the location of the transmitter. Multiple navigation methods are introduced and compared using Monte Carlo Analysis. Despite a noisy signal and a simplified measurement model, the different navigation methods are able to estimate the location of the transmitter with one or more UAVs. Keywords Unmanned aerial vehicle · UAV · Swarm · AggieAir · Potential field navigation · Kalman filter · Radio localization · Fish tracking · Monte Carlo analysis · UAV simulation · Wireless propagation model
1 Introduction In order to preserve our native fish species in the watersheds of the United States, managing, monitoring, and studying fish habitat is very important. Extensive fish habitat studies are conducted to minimize environmental impact before building a dam, while restoring natural habitat in a river, or in ensuring the economic strength of a fisheries market. An important part of these studies is
288
tracking where these fish are migrating, moving and living. This is commonly accomplished by tagging the fish of interest with a radio transmitter then using the transmitter to locate the fish in the future to see where it has moved [4, 5]. Tagged fish are located using a receiver and a directional antenna. As the antenna is pointed in the direction of the fish, the signal sent from the transmitter becomes louder. The operator can hear the beacon through the receiver and navigates toward the fish. While this method has proven to be effective, it can be time consuming and costly. An alternative method for tracking fish utilizes an unmanned aerial vehicle (UAV) developed at Utah State University (USU) called AggieAir [7]. Unlike most UAVs, AggieAir is primarily a remote sensing platform. It carries many types of payload, including a visual camera, a near-infrared (NIR) camera, a thermal infrared camera [8], and air quality sensors. A fish tracking payload has also been recently developed for AggieAir [11]. Using the fish tracking payload, AggieAir is able to fly over a body of water and receive the signal transmitted by the tagged fish. As AggieAir flies over the area of interest, the fish tracker measures the signal strength from its receiver and logs other data, including position, orientation, and speed of the UAV. With this information, the UAV can be directed toward areas of the map to improve the estimate of the tag location. Others have looked into similar scenarios. In a general radio tracking problem, Frew et al. [2] proposed a radio localization technique cast as a distributed estimation problem. Consensus within a group of UAVs was reached with a decentralized method; however, the algorithms used are computationally complex and only considered communication with one neighboring UAV. Korner et al. [10] introduced a method that uses particle filters to track multiple targets with radio tags. Implementation employed only one UAV. While these methods work, a simple, decentralized method is needed for multi-UAV navigation and radio tag localization. A method employing these features could be more practical for use in small UAVs with embedded systems. This paper is an extension of the author’s previous work [6]. This work is focused on simple methods for UAV navigation, multi-UAV
J Intell Robot Syst (2014) 74:287–307
coordination, and radio transmitter location estimation. Simplicity is important for real-world, real-time fish tracking applications. Therefore, swarm-like rules (using potential fields [3]) were established for navigation and multi-UAV coordination/avoidance. Much has been added in this paper. Real data was collected from a UAV and used to develop a propagation model. This propagation model was used to develop a real-world simulation to generate signals for the UAVs. A simplified version of the propagation model and an Extended Kalman Filter were used to estimate the position of the transmitter. Four simple navigation methods to control one or more UAVs are also introduced. These methods generate potential fields based on the estimated location of the transmitter and the covariance matrix of the estimated position. Comparisons are made between the methods by testing them with the simulator and Monte Carlo analysis. The results give clear answers to which methods are effective, and what the expected performance would be in real-world scenarios.
2 Methods In order to properly test different methods to locate tagged fish, a real-world simulation was created. Figure 1 shows a diagram that describes the different pieces of the simulation. Multiple UAVs can be included in the simulation. Each UAV (i) uses wind (|W|, ∠W) and noise (ω, η) to generate ˆ i ). The truth states (Xi ), and estimation states ( X potential field map receives the estimation states from each of the UAVs and updates the map with their current position. If the UAV receives a
Fig. 1 Flight simulation system diagram
J Intell Robot Syst (2014) 74:287–307
289
signal from the fish tag, the UAV estimation states are also used to update the estimated tag position. The truth state of the UAV is used by the tag propagation model to determine whether or not a UAV is close enough to receive a measurement and what that measurement will be ( S˜ i ). When the ˆ T ) is updated with a new estimated tag state ( X measurement, the estimated tag state and the corresponding covariance matrix (PT ) are both used to update the potential field map. The potential field map provides a force (F) which the UAVs will use to control their heading. The UAVs share each other’s position, along with the estimated position and covariance matrix of the tag, indirectly through the potential map. Static objects can also be contained in the potential map like ground based obstacles and a boundary. The tag state estimation error (δ XT ) is used to determine the success of the tag estimation and is calculated by subtracting the tag state estimate with the tag truth state (XT ). Monte Carlo analysis was used with δ XT to determine whether or not a method was successful. The basic idea of Monte Carlo analysis is to run a high-fidelity simulation many times. Because of the noise and disturbances included in the simulation, each trajectory generated the simulation will be different. Monte Carlo analysis is used to show how the trajectories vary and predict expected performance in real-world scenarios. For example, the simulation in Fig. 1 was run 100 times to look at how much δ XT varies for the different methods. One single run by this simulation will be referred to as a flight. Equation 1 shows how an error trajectory δ Xk (t) is calculated for flight number k:
The mean of the error trajectories describes the accuracy that most of the flights should achieve for that given method. The variance describes the best and the worst that can be expected around the mean. The uncertainty of the variance is also considered. Many applications require thousands of runs to reduce the uncertainty of the variance. For this application, 100 flights is enough to compare navigation methods. In addition, the uncertainty of the variance is documented and confidence intervals around the variance curves are calculated using a chi-squared distribution. The subsections below describe the important pieces of the flight simulation. These pieces include the environmental dynamics, the UAV system dynamics, simulating and estimating the tag location, and UAV navigation using potential fields.
ˆ Tk (t). δ Xk (t) = XT (t) − X
b ∠W b˙ ∠W = − + ω∠ W . τ∠W
(1)
After all of the flights are finished, the error trajectories are compiled in (t): (t) = [δ X1 (t), δ X2 (t) · · · δ X N (t)]T .
(2)
Then the mean μ(t) and variance (t)2 of (t) are found for each time (t): μ(t) = E[(t)],
(t)2 = E[(μ(t) − (t))(μ(t) − (t))T ].
2.1 Environmental Dynamics Equation 5 is used to calculate the magnitude and direction of wind. |W|0 and ∠W0 are nominal values of the wind and b |W| and b ∠W are random biases modeled as first order Markov stochastic processes (Eq. 6). |W| = |W|0 + b |W| , ∠ W = ∠ W0 + b ∠ W ,
(5)
˙ 0 = 0, |W| ˙ 0 = 0, ∠W b |W| b˙ |W| = − + ω|W| , τ|W|
The driving noise (ω|W| , ω∠W ) is Gaussian with the following properties. E[ω|W| ] = 0,
E[ω|W| ω|W| ]=
E[ω∠W ] = 0,
E[ω∠W ω∠ W] =
(3)
(4)
(6)
2 2σ|W|
τ|W|
δ(t − t ).
2σ∠2 W δ(t − t ). τ∠W (7)
290
J Intell Robot Syst (2014) 74:287–307
2.2 UAV System Dynamics Figure 2 shows a detailed diagram of the simulated UAV system for U AVi . The force from the potential field map (F) is used by the controller to generate ui (desired rate of change for the heading), which is used along with noise and wind to generate the UAV truth states (Xi ) and GPS measurements (z˜ i ). The GPS measurements are used as measurements in the Kalman Filter to generate the estimated states of the UAV. 2.2.1 Truth Dynamics
Equation 10 shows the equations to calculate other important truth variables: ground speed (vgi ), and heading (φi ). vgi = x˙ 2i + y˙ 2i , φi = tan−1
y˙ i . x˙ i
(10)
Part of the truth model includes generating the GPS measurements that will be used by the navigation. As shown by Eq. 11, the kth measurement is a function of the current position (xi or yi ) and the respective random bias: x˜ i [k] = xi + b xi ,
Equation 8 shows the UAV truth state dynamics. xi and yi describe the position of the UAV. The time derivative of the position is a function of the airspeed vi , the yaw ψi , and the x and y components of the wind (Wx , W y ). The airspeed is assumed to be constant and the change in yaw is equal to the controller output (ui ). x˙ i y˙ i v˙ i ψ˙ i b˙ xi b˙ yi
= = = =
vi cos(ψi ) + Wx , vi sin(ψi ) + W y , 0, ui , b xi =− + ωxi , τx b yi =− + ω yi . τy
E[ω yi ] = 0,
2σx2 δ(t − t ). τx 2σ y2 E[ω yi ωyi ] = δ(t − t ). τy
(11)
A measured ground speed (v˜ gi [k]) and heading (φ˜i [k]) are also generated from the GPS measurements: (x˜ i [k]− x˜ i [k−1])2 +( y˜ i [k]− y˜ i [k−1])2 v˜ gi [k] = , t[k]−t[k−1] φ˜ i [k] = tan−1
y˜ i [k] − y˜ i [k − 1] . x˜ i [k] − x˜ i [k − 1]
(12)
2.2.2 UAV Navigation Estimation
(8)
b xi and b yi are biases in the GPS measurement. These biases are also modeled as first order Markov stochastic processes. E[ωxi ] = 0,
y˜ i [k] = yi + b yi .
E[ωxi ωxi ] =
To control the heading of the UAV according to the potential field map, an estimated position ˆ This is and heading of the UAV is needed (X). achieved through a Kalman filter. ˆ = [xˆ i yˆ i φˆi ]T . X
(13)
To develop the Kalman filter, a design model is used where Xdm = [xdm ydm φidm ]T , i i
(14)
and (9)
dm x˙ dm = v˜ gi cos(φidm ) + ωdm i x , dm y˙ dm = v˜ gi sin(φidm ) + ωdm i y ,
φ˙ idm = ui + ωφdm .
(15)
dm dm ωdm x , ω y , and ωφ represent the process noise of the propagation. Equation 15 can also be represented in the following vector form:
Fig. 2 UAV system diagram
dm ˙ dm = fˆ(Xdm , v˜ gi ˆ dm , X , ui ) + Bw
(16)
J Intell Robot Syst (2014) 74:287–307
where fˆ =
291
where
⎡
⎤ dm v˜ gi cos(φidm ) dm ⎣ v˜ gi sin(φidm ) ⎦ ,
ˆ X ˆ − ), zˆ i [k] = h( ⎡ ⎤ 100 ˆ ∂h Hˆ = = ⎣0 1 0 ⎦ , ˆ ∂X 001
ui ⎤
⎡
100 Bˆ = ⎣0 1 0⎦ , 001 dm
w
=
dm dm T [ωdm x ω y ωφ ] .
(17)
(26) (27) ⎡
⎤ σx2 0 0 Rˆ = E[ν dm (ν dm )T ] = ⎣ 0 σ y2 0 ⎦ . 0 0 σφ2
(28)
Using this design, the filter state is propagated according to the following equations:
2.2.3 Heading Control
˙ˆ ˆ v˜ gi , ui ), X = fˆ(X,
(18)
˙ ˆ Bˆ T , Pˆ = Fˆ Pˆ + Pˆ Fˆ T + Bˆ Q
(19)
All of the UAVs share their location to the potential field map. When U AVi shares its location (xi and yi ), the potential field map finds the gradient of the potential in that area (created from other UAVs, the border, other obstacles, etc) and returns a force Fi (xi , yi ). The x and y components of the force are used to find the desired heading φ f i of the UAV (Eq. 29).
where
⎡ ⎤ 0 0 − sin(φˆ i ) ˆ ∂f Fˆ = = ⎣0 0 cos(φˆ i ) ⎦ , ˆ ∂X 00 0 ⎡ 2 ⎤ σ px 0 0 2 ˆ = E[wdm (wdm )T ] = ⎣ 0 σ py 0 ⎦. Q 2 0 0 σ pψ
(20)
(21)
For the update step of the Kalman filter, Eq. 22 shows how the GPS measurements are related to the states of the design model. dm x˜ dm + ηbdmx , i [k] = xi dm y˜ dm + ηbdmy , i [k] = yi
φ˜ idm [k] = φidm + ηφdm .
(22)
The vector form of Eq. 22 is shown below: ˆ dm ) + ν dm , z˜ i [k] = h(X where
⎡
⎤ dm
xi ⎥ ˆ dm ) = ⎢ h(X ⎣ ydm i ⎦, φidm
(23) ⎡
⎤ dm
ηb x ⎢ ⎥ ν dm = ⎣ηbdmy ⎦ . ηφdm
(24)
Fi (xi , yi ) = Fx ˆi + F y ˆj,
Fy φ f i = tan . Fx
Figure 3 shows a diagram of the closed-loop heading controller. The error (φe = φ f i − φˆ i ) is multiplied by a control gain k in order to give the appropriate amount of control to turn the UAV in the right direction. This control output (ui ) is the turning rate of the UAV and is limited to physical constraints (|ui | ≤ umax ). This physical constraint is applied to the control output before it is given to the UAV. 2.3 Simulating and Estimating the Tag Location ˜ received by the UAVs To simulate the signal ( S) from the fish tag, actual flight data was collected with a UAV to find the propagation model. This data was fit to Eq. 30 where fs is the signal model, r is the distance from the UAV to the tag (xt ,yt ),
Based on this design, the update equations for the Kalman filter are shown below: ˆ −1 , K = P− Hˆ T ( Hˆ P− Hˆ T + R) ˆ − (I − K H) ˆ T + K RK ˆ T, P+ = (I − K H)P ˆ+ =X ˆ − + K(z˜ − zˆ ), X
(25)
(29)
Fig. 3 Closed loop heading controller
292
J Intell Robot Syst (2014) 74:287–307
α is the difference between the polarization angles of the tag and the UAV, and ηs is Gaussian white noise. As Eq. 33 shows, α is symmetric around its own axis and the axis perpendicular to it (0 ≤ α ≤ π/2). S˜ = fs (r, α) + ηs .
(30)
E[ηs ] = 0,
(31)
E[ηsi ηsjT ] = σs2 δ(ij).
r = sqrt(xˆ t − xˆ )2 + ( yˆ t − yˆ )2 ),
(32)
⎧ ⎪ αu − αt + 180 ⎪ ⎪ ⎪ ⎨α − α t u α= ⎪ α − α u t ⎪ ⎪ ⎪ ⎩α − α + 180 t u
.
2.4 Navigation Using Potential Fields 2.4.1 Conventional Potential Fields Potential fields (PFs) are commonly used to dynamically plan the motion of ground [3] and air-based [13] mobile robots. They are used to attract or repel agents to or from objects and other agents. For example, if each agent was surrounded by a repulsive PFs they would all repel each other and avoid collision. Equation 39 shows a conventional definition of an attractive PF [9].
−180 ≤ αt − αu < −90 −90 ≤ αt − αu < 0 0 ≤ αt − αu < 90 90 ≤ αt − αu ≤ 180
(33) An Extended Kalman Filter was used to estimate the position of the tag. Xt is the state vector for the filter. x Xt = t . (34) yt ˙ t = 0); It is assumed that the tag will not move (X therefore, there is no need for a propagation step with this filter. However, a linearized measurement model is needed. Equation 35 shows the measurement model that was used in the Kalman Filter. hˆ is a simplified version of fs and is not a function of α. ˆ Sˆ = h(r).
(35)
With H, the Kalman gain, K, can be computed, and the state and covariance matrices can be updated using Eq. 36: K = P− H T (H P− H T + R)−1 , P+ = (I − KH)P− (I − KH)T + K RK T , ˆ ˆ+ =X ˆ − + K( S˜ − S), X t t
(36)
where, ∂h H= , ∂Xt
(37)
R = E[ηs (ηs )T ] = σs2 .
(38)
U(x) =
1 k(xa − xg )2 . 2
(39)
This PF is applied to the agents through a force defined as the following negative gradient of the PF: F(x) = −∇U(x).
(40)
As xa (position of agent) gets further away from xg (position of a goal), U(x) goes up. The slope of U(x) also increases as U(x) goes up, which increases the force acting on the agent from the PF and attracts the agent toward xg . A repulsive PF would increase as the distance from the object to the agent decreases and would repel the agent. 2.4.2 Fractional Order Potential Fields The method used to generate the PFs for this project is called the Fractional Order Potential Field (FOPF) [1]. The main advantage of the FOPF over other methods is that the FOPF allows the user to design repulsive PFs based on a level of danger. For example, a collision with another UAV would be very dangerous and therefore would get a different PF that produced a higher force than the PF of a less dangerous object like a boundary. A repulsive FOPF is derived starting with the definition of the Coulombian electric field E(r) E(r) =
q . 4π0 r 2
(41)
A single integration of E(r) produces the Coulombian potential from the electric field of a
J Intell Robot Syst (2014) 74:287–307
293
punctual charge (V1 (r)), and a double integration produces the Coulombian PF generated by the uniformly distributed charge along a straight-line segment V2 (r): V1 (r) =
q , 4π0 r
(42)
q ln r V2 (r) = . 4π0
(43)
Both V1 (r) and V2 (r) could be used as repulsive PFs. Figure 4 (for n = 1 and n = 2) shows the normalized potential of V1 (r) and V2 (r) vs. distance and shows that V2 (r) will have a stronger force at a larger distance than V1 (r). For a more dangerous object, V2 (r) would be preferred over V1 (r). It turns out that with each successive integration of E(r), more force is applied at a further distance. Fractional calculus can be used to generate these integrals, as well as fractional order integrals. This conversion using the Weyl fractional integral [12] is given below: Vn (r) = W r E(r) =
q 4π0 (n)
r
∞
(θ − r)n−1 dθ, r2 (44)
Fig. 4 Repulsive Fractional Order Potential Field (1 ≤ n ≤ 5)
where Vn (r) is the nth integral of E(r) (n can be any number greater than 0) and (n) is the Gamma function. 1 n−1 1 (n) = ln dt. (45) t 0 After manipulating Eq. 44, Vn (r) can be written as shown below: Vn (r) =
q (2 − n) ∀n ∈ (0, 2) (2, ∞). 4π0 r 2−n
(46)
To further simplify, Vn (r) can be normalized between 0 and 1 with a maximum and minimum distance (rmax and rmin , respectively) where the normalized potential function U n (r) would be 1 at the minimum distance and 0 at the maximum. Urep (r) =
n−2 Vn (r) − Vn (rmax ) r n−2 − rmax = n−2 . n−2 Vn (rmin ) − Vn (rmax ) rmin − rmax
(47) Figure 4 shows the result of Eq. 47 with 1 ≤ n ≤ 5. When n = 2, the denominator of Eq. 47 goes to zero and produces a singularity; therefore, Eq. 43 is normalized and used in its place. As Fig. 4 shows, an FOPF with a greater n will be more
294
J Intell Robot Syst (2014) 74:287–307
useful in dangerous cases and will apply a greater force sooner than an FOPF with a smaller n. One other advantage of the FOPF over other conventional PFs is that they are normalized and have maximum range. This will keep all PFs equal and unable to affect the agent if it is outside its range. For these reasons FOPFs are also used as attractive PFs (Fig. 5). To find the attractive FOPFs, Eq. 47 is mirrored around the x-axis. U att (r) = −
n−2 r n−2 − rmax n−2 n−2 rmin − rmax
.
(48)
2.4.3 Navigation for Tag Estimation Besides avoiding other UAVs and staying within bounds, the main navigation objective is to fly in areas to improve the estimate of the tag position. Four methods are presented: The Attractive Method, the Repulsive Method, the Offset Repulsive Method, and the Dual-offset Repulsive Method. Figure 6 shows a PF map for the Attractive Navigation method. In this method, a large attractive PF is formed at the estimated position of
Fig. 5 Attractive Fractional Order Potential Field (1 ≤ n ≤ 5)
the tag. This draws the UAVs to the estimated position (the center of the graph in Fig. 6). The Repulsive Method also has a large attractive field centered at the estimated position of the tag. However, it includes a small repulsive PF, also centered at the estimated position of the tag. The attractive field draws the UAVs in, and the repulsive field keeps them from getting close to the estimated position. Instead of going to the estimated position, they circle around it. Figure 7 shows the PF map for the Repulsive Method. The wide band 400 m from the center of the map is where the repulsive field meets the attractive field. This is the area of least potential and the location to which the UAVs will be attracted. The Offset Repulsive Method is identical to the Repulsive Method, except the small repulsive field is centered around a point offset from the estimated position of the tag. Figure 8 shows an example of the Offset Repulsive Method. In this example, the repulsive PF is placed at an offset of 100 m in the positive x direction. This offset creates the moon shaped area displayed in Fig. 8 to which the UAVs are attracted. The axis on which the offset is placed is dependent on the covariance of the estimate. If the offset is currently placed
J Intell Robot Syst (2014) 74:287–307
295
Fig. 6 Potential field map for the Attractive Method
Single Attractive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
200
400
600
800
1000
x(m)
on the x-axis and the UAVs are flying in that area, then the variance of the estimate along the x-axis will decrease. Once it is 10% lower than the variance along the y-axis, the offset will move clockwise to the y-axis and the UAVs will move with it. After a few iterations, the offset and the
Fig. 7 Potential field map for the Repulsive Method
UAVs will have made a full circle around the estimate in a clock-wise direction. The Dual-offset Repulsive Method uses two repulsive PFs. Each one is offset an equal distance from the estimated position of the tag to create the map shown in Fig. 9. In this specific
Single Attractive with Single Repulsive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
x(m)
200
400
600
800
1000
296
J Intell Robot Syst (2014) 74:287–307
Fig. 8 Potential field map for the Offset Repulsive Method
Single Attractive with Offset Repulsive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
200
400
600
800
1000
x(m)
example, the two offsets are placed 200 m away from the estimate on the x-axis. This creates the two low-potential areas shown above and below the repulsive PFs on the y-axis. Like the Offset Repulsive Method, this method also chooses the axis to place the offsets based on the variance of the estimation. If the variance along the x-axis is too high, then the offsets will be placed on the y-axis to attract the UAVs to the x-axis. This method might produce better results with multiple UAVs than the Offset Repulsive Method since
Fig. 9 Potential field map for the Dual-offset Repulsive Method
it encourages the UAVs to spread out and fly in different areas.
3 Results 3.1 Tag Signal Propagation Model To model the signal propagation, the AggieAir fish tracking payload was flown in the flower
Single Attractive with Dual Offset Repulsive PT Map 1000 800 600 400
y(m)
200 0 −200 −400 −600 −800 −1000 −1000 −800
−600
−400 −200
0
x(m)
200
400
600
800
1000
J Intell Robot Syst (2014) 74:287–307
297
the signal strength (approx. 80 dB). This is to be expected in a wireless system, especially while using a relatively simple model. A more complex model may yield lower noise but is not needed for the simulation since the goal is to recreate the signal (even if it is very noisy) and estimate its location with a model that is even more simple. Once the measurement model for the signal was known, hˆ and H from Eqs. 35 and 37 could be specified. As shown in Eq. 50, only the first two terms in r were used for estimation; none of the α terms were used.
Fig. 10 The flower flight path used to collect flight data (footnote 1)
pattern shown in Fig. 101 at 400m above the tag. Before the flight, the position of the tag (xt , yt ) and the polarization angle (αt ) were both noted. During the flight, the payload software logged the signal strength and the position and orientation of the UAV. These data ware used to find r and α according to Eqs. 32 and 33. r and α were then used with the measured signal strength ( S˜ m ) to find the coefficients of fs (Eq. 49). fs (r, α) = β1 + β2r + β3r 2 + β4 r 3 + β5 + β6 α + β7 α 2 + β8 α 3 .
(49)
Figure 11 shows the result of the modeling. The surface in the plot is the model and the points around the surface are the measured signals received during the flight. Figure 12 shows the residuals between the model and the data. The data from this figure show the noise of the signal and are used to generate ηs in Eq. 30 for the simulation. Notice that the 3-sigma value for the noise is half of the maximum possible value for
1 Figure 10 does not show the path of the UAV during the simulations.
ˆ h(r) = β1 + β2r,
(50)
∂ hˆ β2 (xˆ t − xˆ ) β2 ( yˆ t − yˆ ) H= = , . ∂Xt r r
(51)
3.2 Flight Simulation Figure 13 shows the a map of the flight simulation environment. A border surrounds the flight area and has a potential field to keep the UAVs inside. The tag is located in the middle of the map (0,0) and the boundary of the signal is indicated by the black dotted circle around the tag location. The location of the UAVs are indicated on the map by a star symbol. The contours surrounding the UAV are their potential fields, and the blue line behind the UAV is its track. The estimated location of the tag is drawn on the map by the ‘o’ symbol. The large contour lines centered at the estimated tag location are the potential fields of the estimate (Single Attractive Method). The covariance of the tag location estimate is also indicated on the map by the covariance ellipse drawn in black. Initial simulations were run to validate the realworld disturbances added to the simulation (e.g. wind, GPS error, etc). Table 1 contains the values for the simulation parameters that were used for the following results. 3.2.1 Simulated Disturbance and Noise Results Figure 14 shows an example of the position error of the UAV in the x- and y-axis. The simulation did well at creating conditions that a UAV might encounter during flight.
298
J Intell Robot Syst (2014) 74:287–307
Fig. 11 The two-dimensional, fourth-order polynomial model plotted with the test flight data
Actual Signal Data and Model 120 100 Model
Signal
80
Actual Data
60 40 20
0 1000 800 600 r
400 200 0
The wind was also well simulated (Fig. 15). This is even more apparent when looking at the plots of ground speed (Fig. 16) and yaw vs. heading (Fig. 17). As is expected under windy conditions, the ground speed is highly variable even though the airspeed was constant. The difference between the values for yaw and heading (crab
Fig. 12 Residuals of the data vs the signal model
80
100
60
40
0
20
α
angle) also illustrates the success in modeling the wind. Figure 18 shows a plot of the measured signal ˜ vs. the signal strength without noise strength ( S) ( fs ). A noisy, variable signal is expected given the noisy data acquired and the dependence on polarization angle.
Residual error plot. 3−sigma: 40.3372 70 30 60 25
50
20
α
40
30
15
20
10
10
5
0
0
100
200
300
400 r
500
600
700
J Intell Robot Syst (2014) 74:287–307 Fig. 13 Flight simulation map
299 Border Potential Field
1000 800 600
Covariance Ellipse
UAV Location & Potential Field
400 200 0 −200 −400
Estimated Tag Location
Tag Location
−600 −800 −1000
−1000
3.2.2 PF Navigation Method Results The Attractive PF Method is demonstrated in Fig. 13. As expected, the UAV was drawn toward the estimated position of the tag due to the large attractive field centered at the estimated position. Figure 19 shows a flight using the Repulsive Method. With this method, the UAV was encourTable 1 Simulation parameters Symbol
Value
Description
dt GPSdt Sdt 3σx 3σ y 3σφ
0.1 1 4 5 5 2
τx τy |W|0 3σ|W| 3σ∠ W τ|W| τ∠ W umax k xT , yT αT 3σs
10 10 5 5 5 50 50 10 10 (0,0) 45 40.3
Simulation step size (s) GPS Measurement interval (s) Tag interval (s) 3-sigma value for GPS bias (m) 3-sigma value for GPS bias (m) 3-sigma value for heading process noise (degrees) Time constant for GPS bias in x (s) Time constant for GPS bias in y (s) Nominal wind magnitude (m/s) 3-sigma wind magnitude bias (m/s) 3-sigma wind direction bias (degrees) Time constant for wind magnitude (s) Time constant for wind direction (s) Max turning rate for UAVs (deg/s) Gain for heading controller Position of tag (m) Polarization angle of tag (degrees) 3-sigma value for signal strength noise (dB)
−500
0
500
1000
aged to circle around the estimated position; however during a few simulations, it would get turned around by the wind and blocked from covering the entire circle. The Offset Repulsive Method (Fig. 20) was more successful at directing the UAV completely around the estimate. In addition, it yields an interesting variance plot (Fig. 21) that shows “steps” of improvement in each axis when that particular axis is being worked on by the UAV. The Dual-Offset Repulsive Method yielded similar results to the Offset Repulsive Method. Figure 22 shows a flight with two UAVs. The Dual-Offset Repulsive Method was successful at separating the UAVs and directing them toward the two separate areas of low potential. 3.3 Monte Carlo Simulations Table 2 shows the results of the Monte Carlo Simulations. The mean, the 3-sigma values, and the 99 % chi-squared confidence interval for the 3-sigma values are listed for all four methods. Each method was also simulated with one, two and three UAVs. With the highest mean and 3-sigma values, the Attractive method had the worst performance and is not worth looking into further. Focusing on the other three methods and the simulations
300
J Intell Robot Syst (2014) 74:287–307
Fig. 14 Plot of UAV position error in X and Y
Position Error, var x: 2.6454 var y: 1.3555 6 5 4 3
Error (m)
2 1 0 −1 −2 −3 −4
10
20
30
40
50 60 Time(s)
70
80
90
100
70
80
90
100
70
80
90
100
Wind Magnitude (m/s)
Wind Magnitude 10 8 6 4
0
10
20
30
40
50 60 Time(s)
Wind Direction Wind Direction (deg)
Fig. 15 Plot of wind magnitude and direction
0
195 190 185 180 175 170
0
10
20
30
40
50 60 Time(s)
J Intell Robot Syst (2014) 74:287–307
301
Fig. 16 Plot of actual and measured ground speed
Ground Speed 22 Truth Measured
20
Speed(m/s)
18
16
14
12
10
8
Fig. 17 Plot of yaw, heading and wind direction
0
10
20
30
40
50 60 Time(s)
70
80
90
100
Yaw, Heading, Wind Direction Plot 180 160 140
yaw(psi) heading(phi) est heading(phi) Desired phi wind
Degrees
120 100 80 60 40 20 0 −20 5
10
15
20 25 Time(s)
30
35
40
302
J Intell Robot Syst (2014) 74:287–307
Fig. 18 Plot of measured and modeled signal strength
Measured vs Modeled Signal 90 80 70
Signal
60 50 40 30 20 10 0
0
50
with just one UAV, the following observations are noted. 1. The two offset methods have very similar results and have a lower mean value than the Repulsive Method.
Fig. 19 Repulsive method simulation results
100
150 Time(s)
200
250
300
2. The Repulsive Method may have a higher mean than the offset methods, but the upper 3-sigma value is lower. 3. The Repulsive Method is the only one that has a lower 3-sigma value greater than zero.
1000 800 600
wind
400 200 0 −200 −400 −600 −800 −1000
−1000
−500
0
500
1000
J Intell Robot Syst (2014) 74:287–307 Fig. 20 Offset repulsive method simulation results
303 1000 800 600
wind
400 200 0 −200 −400 −600 −800 −1000
−1000
−500
4. The Repulsive has the smallest gap between the mean and the 3-sigma values and the smallest Chi-2 spread.
The conclusion from these observations is that, while the offset methods may achieve higher
Fig. 21 Offset repulsive method tag position variance plot
0
500
1000
accuracy estimates, they are less precise and have higher uncertainty. The Repulsive Method may not be able to achieve errors as low as the Offset Methods; however, is much more reliable. The error trajectory plots of the Repulsive Method and the Offset Method also confirm this conclusion (Figs. 23 and 24).
Tag XY Estimation Variance 900 x y 800
Variance
700
600
500
400
300
200
0
50
100
150
200
250 300 Time(s)
350
400
450
500
304
J Intell Robot Syst (2014) 74:287–307
Fig. 22 Dual-offset repulsive method simulation results (2 UAVs)
1000 800 600
wind
400 200 0 −200 −400 −600 −800 −1000
−1000
The table also shows that, with the Attractive and Repulsive methods, having multiple UAVs improved the estimation performance. The opposite is true with the Offset and Dual-offset methods. For these methods, their accuracy and confidence were reduced as more UAVs were
Table 2 Monte Carlo simulation results for tag error dispersion and 100 flights
All of the results are in meters
−500
0
500
1000
added to the simulation. Even so, the general improvement shown in the table from multiple UAVs is not great. The greatest contribution that multiple UAVs made to the estimate is illustrated in the error trajectory plots. Figure 25 shows the trajectory plot for the Repulsive Method and
Method
# of UAVs
Final mean
Upper 3-sig
Lower 3-sig
Chi-2 spread
Att
1
274
671
0
149
Rep
1
251
343
159
35
Offset
1
160
408
0
93
Dual
1
160
360
0
76
Att
2
249
577
0
123
Rep
2
256
325
187
27
Offset
2
168
411
0
98
Dual
2
186
436
0
101
Att
3
234
572
0
127
Rep
3
255
318
192
27
Offset
3
186
443
0
111
Dual
3
190
471
0
108
J Intell Robot Syst (2014) 74:287–307
305
Fig. 23 MC simulation with 1 UAV using repulsive method
Tag Estimation Error Trajectories 1000 Number of Runs/Errors: 100/0 Final Mean: 251 3−Sig Values, Upper:343 Lower:159 Chi−2 Spread: 34.6246
Estimation Error (m)
800
Mean 3 Sigma 99% Chi2 Bounds
600
400
200
0
−200
−400
0
500
1000
1500
2000
2500
3000
Time (s)
Fig. 24 MC simulation with 1 UAV using offset repulsive method
Tag Estimation Error Trajectories 1000 Number of Runs/Errors: 100/0 Final Mean: 160 3−Sig Values, Upper:408 Lower:−88 Chi−2 Spread: 93.3655
Estimation Error (m)
800
Mean 3 Sigma 99% Chi2 Bounds
600
400
200
0
−200
−400
0
500
1000
1500
Time (s)
2000
2500
3000
306
J Intell Robot Syst (2014) 74:287–307
Fig. 25 MC simulation with 3 UAV using repulsive method
Tag Estimation Error Trajectories 1000 Number of Runs/Errors: 100/20 Final Mean: 255 3−Sig Values, Upper:318 Lower:192 Chi−2 Spread: 26.8739
Estimation Error (m)
800
Mean 3 Sigma 99% Chi2 Bounds
600
400
200
0
−200
−400
0
500
1000
1500
2000
2500
3000
Time (s)
three UAVs. Compared to the Repulsive plot with one UAV (Fig. 23), the plot with three UAVs converges ten times faster than the plot with one UAV. Even though using multiple UAVs did not improve the error in the estimate, they can be used to find the tag in a shorter time.
4 Conclusions A real-world simulation was created to compare the proposed navigation methods in the task of using UAVs to estimate the position of a small transmitter. Comparisons were also made between using single and multiple UAVs. The simulation proved successful at implementing realworld disturbances like GPS bias, a variable ground speed, and a crab angle. The estimation of the tag was also successful despite the noisy signal and the simplified measurement model. The Monte Carlo analysis showed that the Offset and Dual-offset Methods could estimate the tag with the highest accuracy; however the confidence on those estimates was very low. The Repulsive Method was not able to estimate with the same accuracy but was very consistent. The use of multiple UAVs did not improve the estimation accuracy; in some cases, it made the estimates less
accurate. However, multiple UAVs were able to find the tag ten times faster than only one. Future work includes a sensitivity analysis to figure out which error sources are major contributors to the tracking error. With this information, meaningful improvements can be made on the platform to increase the accuracy. Acknowledgements The authors would like to acknowledge Dr. Mac McKee of the Utah Water Research Laboratory. Thank you for your continued support. Gratitude also goes to Miguel Leonardo and Ivan Jimenez for their efforts and skills in developing the fish tracking payload. This work is supported by Utah Water Research Laboratory MLF 2006-2013.
References 1. Chen, Y., Wang, Z., Moore, K.: Optimal spraying control of a diffusion process using mobile actuator networks with fractional potential field based dynamic obstacle avoidance. In: Proc. of the IEEE International Conference on Networking, Sensing and Control, pp. 107–112. IEEE (2006) 2. Frew, E.W., Dixon, C., Argrow, B., Brown, T.: Radio source localization by a cooperating UAV team. In: Infotech@ Aerospace (2005) 3. Ge, S., Cui, Y.: Dynamic motion planning for mobile robots using potential field method. Auton. Robot. 13(3), 207–222 (2002)
J Intell Robot Syst (2014) 74:287–307 4. Gray, R.H., Haynes, J.M.: Spawning migration of adult chinook salmon (oncorhynchus tshawytscha) carrying external and internal radio transmitters. J. Fish. Res. Board Can. 36(9), 1060–1064 (1979) 5. Hillyard, R.W., Keeley, E.R.: Distribution and spawning migrations of fluvial bonneville cutthroat trout in the Bear River, Idaho. Tech. Rep., Idaho State University (2009) 6. Jensen, A., Chen, Y.: Tracking tagged fish with swarming unmanned aerial vehicles using fractional order potential fields and Kalman filtering. In: Proc. of the International Conference on Unmanned Aircraft Systems (ICUAS) (2013) 7. Jensen, A.M., Hardy, T., McKee, M., Chen, Y.Q.: Using a multispectral autonomous unmanned aerial remote sensing platform (AggieAir) for Riparian and Wetland applications. In: Proc. IEEE Int. Geoscience and Remote Sensing Symp. (IGARSS) (2011) 8. Jensen, A.M., McKee, M., Chen, Y.: Calibrating thermal imagery from an unmanned aerial system— AggieAir. In: Proc. IEEE Int. Geoscience and Remote Sensing Symp. (IGARSS) (2013)
307 9. Khatib, O.: Real-time obstacle avoidance for manipulators and mobile robots. In: Proc. of the IEEE International Conference on Robotics and Automation, vol. 2, pp. 500–505. Institute of Electrical and Electronics Engineers (1985) 10. Korner, F., Speck, R., Goktogan, A.H., Sukkarieh, S.: Autonomous airborne wildlife tracking using radio signal strength. In: Proc. of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 107–112. IEEE (2010) 11. Leonardo, M., Jensen, A., Coopmans, C., McKee, M., Chen, Y.: A miniature wildlife tracking UAV payload system using acoustic biotelemetry. In: Proc. of the ASME International Design Engineering Technical Conferences & Computers and Information in Engineering Conference (2013) 12. Miller, K.: The Weyl fractional calculus. In: Ross, B. (ed.) Fractional Calculus and its Applications. Lecture Notes in Mathematics, vol. 457, pp. 80–89. Springer Berlin Heidelberg (1975) 13. Paul, T., Krogstad, T.R., Gravdahl, J.T.: Modelling of UAV formation flight using 3D potential field. Simul. Model. Pract. Theory 16(9), 1453–1462 (2008)
