Positioning and tracking of
target drone
Anna Hanström and Jet Verheij LiTH-ISY-EX--21/5362--SE
Supervisor: Senior Scientist Rolf Gustavsson
Swedish Defence Research Agency (FOI) Research Engineer Andreas Höggren
Swedish Defence Research Agency (FOI) PhD student Özgecan Özdogan
isy, Linköping university
Examiner: Docent Danyo Danev
isy, Linköping university
Division of Communication Systems Department of Electrical Engineering
Linköping University SE-581 83 Linköping, Sweden
This master thesis studied methods for tracking and localising a moving target from an autonomous seeker drone. Feasible methods for automatic control of the seeker drone and different antenna configurations were explored as well. Two different tracking filters and two different controllers were tested for this pur-pose. The algorithm was developed in Python and MATLAB. The evaluation of the filters and controllers was done both theoretically with simulations but also practically with flight tests. Performance and robustness were measured by exam-ining the estimated target position and the smoothness of the seeker path. Both filters performed satisfactorily, the same conclusion could be made for the auto-matic controllers as well.
Regardless of the sufficient results, for future work there are several aspects which can be improved. The communication with the drone’s motors, the noise models and one of the automatic controllers are all examples of areas which can be improved further.
Firstly, we would like to thank our supervisors at FOI, Rolf Gustavsson and An-dreas Höggren. Without your constant support and great knowledge this project would not have been what it is today. Additionally, we would also like to thank our supervisor and examiner at Linköping University, Özgecan Özdogan and Danyo Danev. We would like to thank Simon Ahlberg, Johannes Lindblom and Patrik Hedström at FOI, your helpful insights always managed to keep us on the right track. Lastly, we are very thankful for the help we have received from Joakim Rydell, Alfred Patriksson and the rest of the people at FOI who have con-tributed to this master thesis in one way or another.
Linköping, February, 2021 Jet Verheij and Anna Hanström
Notation ix 1 Introduction 1 1.1 Motivation . . . 1 1.2 Purpose . . . 1 1.3 Problem statement . . . 1 1.4 Limitations . . . 2 2 Theory 3 2.1 Positioning methods . . . 3
2.1.1 Direction finding method and positioning method - what is the difference? . . . 3
2.1.2 Fundamental ambiguity . . . 4
2.1.3 Received Signal Strength Indication (RSSI) . . . 5
2.1.4 Time difference of arrival (TDOA) . . . 8
2.1.5 Interferometry . . . 10 2.2 Tracking . . . 12 2.2.1 Model formulation . . . 12 2.2.2 Error sources . . . 13 2.2.3 Filters . . . 14 2.2.4 Resampling methods . . . 16
2.2.5 UAV flight controller . . . 19
2.2.6 Optimality criterion . . . 20 2.3 Noise distribution . . . 21 2.3.1 Student-t distribution . . . 21 2.3.2 Rice distribution . . . 22 3 Method 25 3.1 System overview . . . 25 3.2 Setup . . . 27
3.2.1 Global coordinate system setup . . . 27
3.2.2 Global and local coordinate system setup . . . 28
3.2.3 Conversion of GPS data to Cartesian coordinates . . . 28 vii
3.3 Tracking . . . 28
3.3.1 Kalman Filter . . . 29
3.3.2 Particle Filter . . . 29
3.4 Flight control . . . 30
3.4.1 MDP controller . . . 32
3.4.2 Relative bearing controller . . . 32
3.4.3 Collision handling . . . 32
3.5 Test framework setup . . . 33
3.5.1 Test model for signal strength . . . 33
3.5.2 Test model for bearing . . . 37
3.6 Practical tests . . . 38
3.6.1 Simulation of bearing and signal strength . . . 38
3.7 Hardware setup and software programs . . . 38
4 Results 41 4.1 Coordinate system setups . . . 41
4.2 Antenna configuration and positioning method . . . 41
4.3 Tracking . . . 42
4.3.1 Kalman filter . . . 42
4.3.2 Particle filter . . . 45
4.4 Flight control . . . 48
4.5 Noise sensitivity analysis . . . 49
4.5.1 Large bearing noise . . . 50
4.5.2 Large signal strength noise . . . 54
4.5.3 Medium large bearing noise and signal strength noise . . . 58
4.6 Practical tests . . . 61 4.6.1 Initial tests . . . 61 4.6.2 Qualitative performance . . . 64 4.6.3 Flight tests . . . 67 5 Discussion 75 5.1 Signal models . . . 75
5.2 Coordinate system setup . . . 75
5.3 Filters . . . 76
5.4 Controllers . . . 77
5.5 Practical flight test . . . 77
6 Conclusion 79 6.1 Problem statement . . . 79 6.1.1 Positioning method . . . 79 6.1.2 Tracking method . . . 80 6.1.3 Antenna configuration . . . 80 6.1.4 Flight control . . . 80 6.2 Future work . . . 80 Bibliography 81
Abbreviations
Abbreviation Meaning
3d 3 Dimensional
cdf Cumulative Distribution Function
crlb Cramér Rao Lower Bound
foi Swedish Defence Research Agency
gps Global Positioning System
los Line of Sight
mdp Markov Decision Process
pf Particle Filter
rmse Root Mean Square Error
rssi Received Signal Strength Indicator
rt90 Rikets Koordinatsystem 1990
snr Signal to Noise Ratio
tdoa Time Difference of Arrival
uav Unmanned Aerial Vehicle
vhf Very High Frequency
wgs84 World Geodetic System 1984
1
Introduction
In this chapter the motivation and purpose of this master thesis work are pre-sented. The problem statement and the delimitations connected to them are de-fined as well.
1.1
Motivation
In the last couple of years there has been an increased amount of incidents where unmanned aerial vehicles (UAVs) have violated or intruded unauthorised areas. These violations have been especially problematic for airports and public build-ings. Due to these incidents it is of great interest to investigate what methods there are to detect UAVs and localise them in an efficient and rapid way [31].
1.2
Purpose
The focus of this master thesis was to investigate several methods for Unmanned Aerial Vehicle (UAV) detection and localisation using radio frequency signals. The detection and localisation aimed to be implemented on a seeker drone which in a finalised stage would operate and follow the intruding drone completely au-tonomously.
1.3
Problem statement
• What positioning method is most suitable for this task: received signal strength, interferometry or Time Difference of Arrival (TDOA)?
• What method is most suitable for tracking the target drone? 1
• What antenna configuration is the most suitable for this task?
• How should the automatic controlling of the seeker be implemented in or-der to make direction changes smooth?
1.4
Limitations
The task was delimited to a two dimensional case where the height of the seeker drone and the target drone was omitted. Also, the drones were assumed to be in Line of Sight (LOS) from each other, meaning that there would be no obstacles be-tween them. These limitations were made to simplify the tracking and automatic controlling of the seeker drone since less parameters had to be taken into account. When performing practical tests the seeker drone did not have any antenna hard-ware for receiving the signals from the target drone. The focus was on testing the tracking algorithms and automatic controllers. Instead, the Global Positioning System (GPS) on the drones provided the position and from that the bearing and signal strength were simulated, together with noise.
The power of the outgoing signal strength from the target drone was assumed to be known and constant. The reason to have a known power was that it sim-plified the calculations of the noise when the signal strength was simulated. Yet again, this was done to keep the focus on the tracking and automatic controllers. A more elaborated noise model could have been used, but that would have been on the expense of a less studied tracking and/or controller algorithm.
2
Theory
This chapter provides relevant theory for the topics of this master thesis. Section 2.1 first explains the fundamental ambiguity and thereafter consists of the posi-tioning methods: received signal strength, interferometry and TDOA. Section 2.2 explains different methods for tracking a moving object and mentions the filters used in this project. Section 2.3 provides theory regarding the simulated noise distributions that were used for the signal strength and the bearing.
2.1
Positioning methods
This section presents relevant theory regarding positioning methods.
2.1.1
Direction finding method and positioning method - what is
the difference?
A direction finding method gives the direction to the object. A direction can be measured in bearing, the angle in the horizontal plane between true north and the measured direction [27]. A positioning method estimates the exact position and can be measured in multiple ways, for example by the x, y, z position, the bearing together with a distance or with GPS coordinates.
In many cases the method of finding the direction is the first step of the posi-tioning method [3]. For the tracking algorithm it is sufficient to only know the bearing. In this project the distance between the drones was used to avoid colli-sion. Still the tracking algorithms distance estimation properties were examined for both long and short distances, this was to answer the first and second problem statements.
2.1.2
Fundamental ambiguity
The goal was to calculate the exact position of the target drone. For this the bearing α of the target drone compared to the seeker drone had to be estimated, see Figure 2.1. α[n] Target at time step [n] Seeker at time step [n] Global x axis Global y axis
Figure 2.1: A depiction of how the bearing α relates to the seeker and the
target drone.
Here the bearing angle α denotes the angle to the target with respect to the seeker heading.
The first step in that estimation was to receive the signal from the target drone. This could be done with an antenna setup placed on the seeker drone. One gen-eral problem, independent of chosen method, is ambiguity of direction. The am-biguity can be reduced by using different kinds of antennas, different amounts of antennas or by limitations on the problem. One limitation is already mentioned in the section above, Section 2.1.1,"A direction can be measured in bearing, the an-gle in the horizontal plane between true north and the measured direction". Here the horizontal plane was a limitation. In this section the fundamental ambiguity problem will be discussed, in the following sections the three different methods of interest will be explained together with methods to solve the ambiguity. In the simplified case of estimating the direction of a planar wave it is not enough to only recognise the field with a simple antenna in one single point. The
infor-mation must be compared to at least one other point to say something about the direction. When comparing the field in two points it is possible to estimate the angle of direction relative the line connecting the points. In the 3 dimensional (3D) case all the possibilities of direction generates something that could be de-scribed as a cone, this is shown in Figure 2.2 [3]. To remove some ambiguity in a 3D world at least three points are needed, reducing the ambiguity to two, one direction above and one direction below the plane spanned by the three points [3].
Figure 2.2:When measuring a planar wave in two points the different
possi-ble directions create a rotational concentric ambiguity in the 3D case. In the 2D case there will be an ambiguity of two directions.
A commonly used simplification is to do the estimations in the horizontal plane but keeping two simple antennas. This reduces the ambiguity to two directions, one possibility on each side of the line.
2.1.3
Received Signal Strength Indication (RSSI)
The RSSI method studies the signal strength in different positions or in different directions and uses this information to estimate the bearing. The mathematical setup depends on the used method, in this section some general methods will briefly be discussed. Additionally, two different types of antennas will be shortly explained.
Antennas for RSSI
For the purpose of transmitting or receiving a signal there are different antennas which can be used. For the RSSI cases which are described in this section, only directional antennas were used. A directional antenna transmits and receives sig-nals in a specific direction, the direction is dictated by an angle γ which denotes the beam width. The area covered by a directional antenna resembles the shape of a cone. However, there is another type of antenna worth mentioning as well
-an omnidirectional -antenna. An omnidirectional -antenna receives or tr-ansmits in a circular pattern, similar to a bagel, thus covering a wider area than a directional antenna [30].
γ
Figure 2.3: An example of the area which a signal propagates within for a
directional (to the left) and an omnidirectional antenna (to the right).
The RSSI method
The RSSI method builds on the attenuation of the signal, it is a bit like our hear-ing, we turn our head to the direction where the signal is louder. The attenua-tion of the signal depends on the radio propagaattenua-tion channel, this is a combina-tion of the signal frequency and the propagacombina-tion environment. Under ideal cir-cumstances with no obstacles or refraction the attenuation is lower than if there would be obstacles such as trees and bushes. This assumes a free space propaga-tion of the signal, where the signal only takes one path. For ideal condipropaga-tions the signal level reduction loss is described as follows
Ideal condition: 20 · log10(dinitial/df inal) (2.1)
implying that when doubling the distance the received signal level is reduced by 6 dB. Egli [11] proposed an empiric channel model for VHF transmission at ground level stating that the signal level was reduced according to the following equation
Egli: 40 · log10(dinitial/df inal) (2.2)
which implies that by doubling the distance the received signal level would be reduced by 12 dB. Several meters up the condition was closer to ideal, but if the drone was flying on a low altitude the received signal would be much lower. This equation assumes a multi path propagation of the signal, which in practice means that the propagated signal takes more than one path. The cause of this is the environmental conditions which give rise to reflections. The multiple paths introduce distortion which become even more serious when the transmission is in a situation of movement. To redeem the impact of the multi path distortion one can use directional antennas, which has been proved to reduce it greatly [11].
1
2
Figure 2.4:The triangulation method. Two rotations generate two bearings.
Calculating the position where the two bearings cross each other gives the position of the radio source.
RSSI with directional antennas
InMethods for locating signal jammers with a UAV [16] Höggren and Lindmark
used one single directional antenna mounted underneath a UAV. The directional antenna received a higher signal strength when pointing at the jammer. By rotat-ing the UAV a full cirle two times on different places the position of the jammer was identified. This is called the triangulation method, see Figure 2.4. After the first rotation and bearing estimation the drone flew to another position perpen-dicular to the bearing and rotated another time. The two bearings were used to estimate the position. This method generated a position in a fairly straight-forward way, a drawback was that the rotations took time and the calculation depended on an immobile target. If the target drone moved this would not be a useful method for the situation. Also, the method was sensitive to different weather conditions, the full rotation took some time to carry out and the wind could cause it to tilt, rotate at uneven speed and gain or lose altitude.
InHunting Drones with Other Drones: Tracking a moving Radio Target [22] Dressel
and Kochendfer used two directional antennas placed on the opposite sides of the seeker drone, one pointing along the front of the drone, the other one point-ing to the back. Like before the directional antennas received a higher signal strength when pointing at the radio source. By having two directional antennas it was possible to say if the target drone was in front or behind the seeker drone, not if the target drone was to right or left hand side of the seeker drone. This is the fundamental ambiguity and Dressel and Kochendfer used a tracking method described in Section 2.2.6 to solve the problem.
RSSI with multiple antennas
InElectronic warfare using a software defined platform [28] six VHF antennas were
used in different positions together with one signal transmitter. Combining the receivers two and two it was possible to find an oval for each pair, combining the three ovals from the pairs they all intersected in one point, namely the position of the transmitter.
If a pair of receivers got a difference in signal level to, for example, 6 dB and the conditions were ideal, the distance between the emitter and receiver differ-entiated with a factor of 2. This meant that for receiver A the transmitter was twice as far away as from receiver B. This is shown in Figure 2.5 where the blue points indicate the receivers and the red circle indicates possible positions for the transmitter.
A B
Figure 2.5:An example of how the ovals are created, the receiver A got 6 dB
lower signal level than receiver B, implying the receiver is twice as far away from A than from B. The red oval depicts all possible positions where the transmitter possibly can be.
2.1.4
Time difference of arrival (TDOA)
In general there are two different techniques for using TDOA, one of them is where the subtraction between two different measured time arrivals is calculated. The other method uses the cross-correlation where the measured arrival time for one sensor is correlated with the measured arrival time for another sensor. The most commonly used method is the one which makes use of the cross-correlation [1]. Both of these methods are described shortly below.
Subtraction of different measured time arrivals
In a simple case, with no noise and two half wave dipole antennas with distance d between them, the angle of incidence is defined such that
α = −arcsinτ21
d/c. (2.3)
Here τ21 = t2−t1 is the difference in time between the second antenna and
the first antenna receiving the same signal. c is the speed of light [3]. Figure 2.6 shows the formation and how α should be interpreted.
d d sin
α α
α
Figure 2.6: The two antennas are placed with the distance d and α is the
angle of the direction of the incoming planar wave.
Apart from the fundamental ambiguity the estimation method is unambiguous, however, the method can be hard to implement if the distance between the
an-tennas is small. A small distance between the anan-tennas generates a very short τ21
which requires a very accurate clock [3].
Equation 2.3 originates from that t1and t2are known, this is not always the case.
If the incoming wave would be single pulses, and the time between the pulses is long enough (or d is large enough) to not mix the pulses up, it is pretty straight forward to calculate the time differences. If not, which occurs most of the time, multiple samples of the received signal must be used, compared and averaged to estimate τ21.
Cross-correlation
A cross-correlation looks at the similarity of two different signals. One signal is denoted as a reference signal whereas the other is the same signal but including a time displacement as well. The signals are measured by two different sensors, where the reference signal originates from a reference sensor. The displaced sig-nal originates from another sensor, located at a known distance d from the
ref-erence. The measurement at the reference sensor x1, containing the signal s can
then be denoted as
x1(t) = A1· s(t) + n(t). (2.4)
The corresponding measurement x2at the other sensor, which contains the time
displacement τ, can then be described as follows
x2(t) = A2· s(t − τ) + n(t). (2.5)
For these equations n(t) denotes measurement noise, A1 and A2 are signal
am-plitudes. The objective is to find the time displacement τ, which corresponds to the relative time difference of arrival between the two sensors. The displace-ment is obtained by correlating the two signals for different displacedisplace-ments. The
correlation peaks at the location for the correct displacement, thus the estimated time difference τ. The cross-correlation can be summarised into the following equation ˆ Rx1,x2(m) = 1 N N −1 X n=0 s(m)s(n + m) (2.6)
where N represents number of samples [1]. The angle α can then be calculated according to Equation 2.3.
2.1.5
Interferometry
When the distances d between the antennas must be short, and the TDOA there-fore can be hard to implement, interferometry can be a good alternative. A draw-back of interferometry is the troublesome ambiguity situation if d is greater than half the wavelength, λ/2, of the incoming signal. This is in addition to the fun-damental ambiguity. In Hunting Drones with Other Drones [22] they used a frequency of around 900 MHz, and therefore a wavelength of approximately 33 centimetres and half a wavelength of 16.5 centimetres. Depending on the size of the drone and used frequency it is possible to have a longer distance between the antennas than half the wavelength and the ambiguity requires therefore further explanation.
Interferometry uses the phase difference between the received signal in only one
sample. The true phase difference, Ψ21, is defined as follows
Ψ21= 2πdsinα
λ (2.7)
given that two antennas are used. Ψ21 can both be positive, negative and
several multiples of 2π if the quota d/λ ≥ 0.5 . The phase difference, φ21, is easier
to estimate by the following equation
φ21= Ψ21+ n · 2π (2.8)
where n is an integer that gives φ21a value inside the interval −π < φ21≤π.
Combining Equation 2.7 and Equation 2.8 gives Equation 2.9 for the direction of the incoming wave can be defined as
α = arcsin φ21 2π −n d/λ . (2.9)
Keeping the distance between the antennas smaller than half a wavelength is a useful method to keep the ambiguity minimal, but a drawback is that the accu-racy of the phase detector decreases with decreasing d. To avoid ambiguity but keeping a high accuracy it is possible to add an extra antenna, the distance be-tween the two first antennas will be smaller than half the wavelength and the distance between the last and the first antennas will be longer. This is shown in Figure 2.7.
Still there is the problem of the fundamental ambiguity of not knowing if the correct bearing is the estimated angle or the reflection of it over the antenna ar-ray (assuming a 2D environment). Like before this problem can also be solved by adding more antennas. If unambiguity is of higher interest than accuracy (or if the antennas are placed with a short distance) it is enough to add just one an-tenna as shown in Figure 2.8a. Here the third anan-tenna must not be placed on the same line as antenna one and two. The thick line indicates the base pairs, antenna one together with antenna two and antenna two together with antenna three. Together these pairs estimate four directions, two of them in the same direction indicating that that is the right angle.
d d1 sin α α α d1
Figure 2.7: The line-up uses three antennas to make sure the ambiguity is
minimal but keeping the accuracy high.
In the case of a high demand on both accuracy and unambiguity it is possible to expand to five antennas such as shown in Figure 2.8b. For symmetrical reasons it is common to create a 90 degree angle between the line-ups.
1 2
3
(a)When using three anten-nas the fundamental ambi-guity is removed.
(b) Using five antennas re-moves the fundamental am-biguity and generates the possibility to increase the ac-curacy.
Figure 2.8:Two different antenna line-ups.
2.2
Tracking
This section presents the required theory for tracking.
2.2.1
Model formulation
The first step which is needed for tracking is to identify a model which represents the problem. For this kind of task a state space model is a common approach to formulate a model. A general formulation of the state space model can be described as
sk+1= f (sk, uk, wk) (2.10)
yk= h(sk, uk, ek) (2.11)
where sk denotes the state at time instance k, uk is a known control input, yk
denotes the measurement at time instance k. wkand ek are process and
measure-ment noise vectors. In some cases the control input, uk, is omitted depending
on the setup. sk holds data from the dynamic model. f denotes the state matrix
and h denotes the observation matrix [15]. Depending on the hardware setup the
state sk could include the bearing and signal strength. This is what was assumed
for this project, the state formulation could then be summarised to sk= [βk, vk]T
where βk denoted the relative bearing to the target with respect to the seeker
drone’s direction at time instance k [21]. vk represented the signal strength from
the target for time step k. Other possible additions to the state is the distance between the seeker and the target at time instance k and the seeker’s heading for each time instance. Furthermore, the relative velocity of the seeker with respect to the target at time each instance could also be of interest. This representation makes sure that the movements of the target and data related to it are being tracked properly [13].
2.2.2
Error sources
One important error source is the noise, since hardware equipment and process noise may arise it is crucial to know how to handle them. A process noise in this
context is an uncertainty which drives the state, sk+1, forward.
Firstly it is necessary to classify the noise, depending on if it’s white or coloured, Gaussian or non-Gaussian. The difference between a case where white noise is included compared to a case with coloured noise can be vastly different [25].
GPS offset
Introduction to GPS: The Global Positioning System, Second Edition [12] discussed clocks and the error generated by them. Each GPS satellite contains out of four atomic clocks, two cesium and two rubidium. Only one of them is used, the other three are for backup. The clocks are accurate, but not perfect. The satellite clock error is somewhere between 8.64 nanoseconds and 17.28 nanoseconds per day. This generates a corresponding range error of 2.59 meter to 5.18 meter.
InA sensitivity analysis of extended Kalman filter for GPS position estimation with and without clock offset [10] a test was performed where the different mean values of x-position error, y-position and z-position error were calculated for 22 hours. The mean for the x-position error was 1.28 meter, the mean for the y-position error was 2.34 meter and for the z-position error was 2.72 meter. Here it was also shown that knowledge of the clock error could minimise the error further. The GPS offset affected this project when the algorithm was tried in practise. Due to the fact that the outgoing signal from the target was simulated from the GPS coordinates, a few meters off could generate problems mainly when the drones were close to each other.
Path loss factor
A sensor setup in practice rarely or never performs perfectly, the received mea-surements therefore include noise and error sources, such as the Doppler effect. Subsequently it is convenient to model the RSSI with these error sources in mind to obtain the correct measurement [25]. The path loss factor is an important fac-tor when modelling the propagation of the received signal strength. Moreover it is necessary when the distance to a target is estimated. The path loss represents the attenuation between the transmitter and the receiver. The factor is related to environmental conditions as well as the hardware setup and subsequently, de-termines the intensity of the received signal strength. If the hardware setup and the environmental conditions are known the path loss factor can be predicted [2]. As the seeker drone moves along a certain path the environment changes with it, this constant change affects the path loss factor as well. Therefore it is of great importance to model the path loss according to these changes and then apply it
to the measured data [25]. A generalised equation for the path loss, Lp, can be formulated as follows which is the Friis equation formulated in logarithmic units
Lp[dB] = Pt[dBm] − Pr[dBm] + Gt[dB] + Gr[dB] (2.12)
where Ptand Prdenote the transmitter and receiver powers. Grand Gtdenote the
receiver and transmitter antenna gains [2]. As already mentioned the path loss is of importance when measuring the signal strength with antennas. In this project
however, no antennas were used and consequently Gr and Gtwere omitted from
the project scope.
2.2.3
Filters
In order to make a prediction for each time instance a filter is needed. Two dif-ferent filters were implemented in this project, the theory regarding them is pre-sented in this section.
Kalman filter
The Kalman filter is a sequential recursive minimum mean square error estima-tor that allows the estimation of s[n], based on the data {x[0], x[1], ..., x[n]} as n increases. This filter estimates ˆs[n], which is based on the previous estimation ˆs[n − 1]. The filter’s state equation and the filter’s observation equation are as follows
s[n] = as[n − 1] + u[n] (2.13)
x[n] = s[n] + w[n] (2.14)
where u[n] and w[n] are zero mean Gaussian noise with independent samples and E(u2[n]) = σu2 and E(w2[n]) = σw2. a is the control input, if there is any pre-known information about the situation a can be used for this. For example, if the Kalman filter is used for estimation of the height of a drone when it is landing and it is known that for each sample the drone should drop approximately 25%, a = 0.75 would be reasonable. Furthermore, it is assumed that s[−1] = N (µs, σs2). A summation of the Kalman filter equations are shown below [19]
Prediction:
ˆs[n|n − 1] = a ˆs[n − 1|n − 1] (2.15)
Minimum Prediction Mean Square Error:
Kalman Gain: K[n] = M[n|n − 1] σw2+ M[n|n − 1] (2.17) Correction: ˆs[n|n] = ˆs[n|n − 1] + K[n](x[n] − ˆs[n|n − 1])) (2.18)
Minimum Mean Square Error Matrix;
M[n|n] = (1 − K[n])M[n|n − 1] (2.19)
where the indexes are as follows [time index | measurement index]. M[n] is the
minimum mean square error when s[n] is estimated based on the previous data and K[n] denotes the Kalman gain.
The first time the filter runs initial values for ˆs[n − 1|n − 1] and M[n − 1|n − 1] are chosen.
Particle filter (PF)
One type of filter which has been widely used for localisation and prediction of moving platforms is the particle filter, PF. The PF provides an approximation of
the state skbased on a set of samples, so called particles which all have individual
weights, wk [15]. In this project the state sk denoted the relative bearing and the
signal strength. A general outline of the PF is stated below [4].
Algorithm 1Basic outline of a particle filter
InputParticle weights wk−1, model state sk−1, measurements yk
OutputNew model state s∗k, particle weights wk
1: for all particles i do
2: Draw a particle i from sik
3: Assign particle i a weight wik
4: end for
5: Calculate total weight t =P
iwik
6: Normalise weights wki according to t−1wik
7: Calculate ˆNef f =P wik2
8: if ˆNef f less than threshold then
9: Resample
10: end if
11: Estimate new model state s∗k
The assigned weight wki for particle i is defined as follows wk|ki ∝wi
k−1|k−1p(yk|sk−1i ) (2.20)
where yk is the measurement, wk−1|k−1i is particle i’s weight from the previous
time step k − 1. p(yk|sk−1i ) is the probability density function. Index i denotes
the particle i and k denotes the time and measurement index [14]. The drawing
of particles from sik was generated from the state model and was calculated as
follows
sik= a · sik−1−b · ˙βk (2.21)
where ˙βk denoted the measured velocity for the relative bearing β. a and
b were constants specified for the model. The measurement equation was then calculated accordingly
yk= sik. (2.22)
The calculation of ski followed the same pattern for the signal strength. The
equation describing the drawing was outlined as follows
sik= sk−1i − ˙vk (2.23)
where ˙vkdenoted the measured velocity for the signal strength. The
measure-ment ykwas then calculated such that
yk= sik. (2.24)
2.2.4
Resampling methods
There are different methods for resampling particles of a Particle filter, however the main objective is to eliminate particles with low probability and maintain those with high probability. Depending on the application some resampling methods work better than others. It is important that the resampling keeps the balance between maintaining diversity and at the same time keeping the compu-tational complexity down. There are in general four basic resampling strategies; multinomial resampling, stratified resampling, systematic resampling and resid-ual resampling [18]. The general outline of these methods is presented below. All of the resampling strategies presented in this section are modifications of a weighted approximate density. The modifications generate an unweighted den-sity. The reason to make it unweighted is due to the elimination of particles with low probability. As particles of low importance are removed, particles of high importance are multiplied. Thus, the indices in a vector where low probability particles were held are replaced with particles of high probability. This
state-ment can be summarised into the following equation, where pN is the weighted
pN(x) = N X
i=1
wiδ(x − xi) (2.25)
and is replaced by ˆpN, which can be described as
ˆ pN(x) = N X k=1 1 Nδ(x − x ∗ k) = N X i=1 ni Nδ(x − xi), (2.26)
here ni is the number of copies of the high probability particle xi [17].
Multinomial resampling
The multinomial resampling method is one of the most basic resampling meth-ods. Simply put, it consists of drawing new indices from [1, N ] where N is a non-random integer. The drawing acts independently of the existing
approxi-mate weighted density pN. A more hands on implementation of this strategy is
presented below [8].
1. Draw an independent, random uniform number uk∼Ui1≤i≤nfrom an
inter-val (0, 1]
2. Assign index i a value according to Qi−1< uk≤Qi where Qi=Pis=1ws 3. Select xk∗ according to x∗k= x(F−1(uk)) = xi
Here F−1(u
k) denotes the inverse of the cumulative distribution function based
on the normalised particle weights [17]. x∗k is the new resampled particle state.
Item 3 in the list above is utilised for all subsequent resampling methods pre-sented below. The main difference between these resampling techniques is how
the uniform number ukis selected. The resulting drawing of ukand how it differs
depending on the resampling method is illustrated below in Figure 2.9.
Stratified resampling
Stratified resampling consists of generating N subsets in the interval [0, 1). ukis
then drawn independently in each of these subsets [8]. The method for selecting
the elements in uk is summarised in the following equation
uk=(k − 1) + ˜uk
N , where ˜uk∼U [0, 1) (2.27)
where N denotes amount of ordered random numbers. The elements in uk
will all lie separately within their respective subset, thus one subset will only be
occupied by one element. Furthermore, all elements in uk has a value inside of
the range [0, 1). Thereafter x∗k is selected according to item 3 in the list from the
paragraph about the multinomial resampling, which is presented above. Conse-quently, xk∗ is the new resampled particle state.
Systematic resampling
As for systematic resampling it follows the same pattern of generating N subsets as for stratified resampling. Namely, it generates subsets in the interval (0, 1] and
draws ukindependently according to Equation 2.27. However, systematic
resam-pling deviates from the stratified resamresam-pling by the means of that the generated
elements in ukare not randomly ordered. The difference between systematic and
stratified resampling is depicted in Figure 2.9, which is presented below.
When the procedure of drawing uk is finalised, x
∗
k is then selected according to
item 3 in the list from the paragraph about the multinomial resampling, which is
presented above. x∗k is now the new resampled particle state.
Residual resampling
Residual resampling differs slightly from the previously mentioned resampling techniques. This resampling strategy can be combined with any of the other pre-sented resampling techniques above. In residual resampling, the first step is to
determine how many copies, n0i, of a particle xi that should be retained. This is
done by using its weight, wi, such that
n0i = bN wic (2.28)
where N is the total number of particles. The brackets in Equation 2.28 take
the integer part of N wi. Thence, a portion of these n
0
i particles are resampled.
How large the portion is is determined by m, which is defined as
m = N −Xn0i. (2.29)
The resampling is then executed such that a new set of copies, n00i, is generated
where the probability to select xi equals the new weight w
0
i = N wi−n
0
i. In other
words, w0i is assigned the decimal part of N wi. The selection of xi can be
per-formed using any of the presented resampling methods above. More plainly
speaking, wi0 replaces ws in item 2 from the list from the paragraph about the
multinomial resampling. Consequently xi replaces x
∗
k item 3. uk is selected
ac-cording to the chosen resampling technique.
A more illustrative demonstration of how these different resampling methods work is presented below in Figure 2.9. The dots in the figure illustrate how the
uniform numbers in ukare selected depending on the resampling method. From
the figure it is evident that the stratified and systematic resampling have similar
subsets for selecting uk. However, it is also clear that the stratified method picks
uk randomly within each subset. This is not the case for the systematic
resam-pling. The multinomial resampling draws ukrandomly in the interval [0, 1). This
is also evident in the picture since the subsets seem to have no effect on the draw-ing.
Systematic resampling was implemented for this project in order to eliminate the particles with low probability. The reason why systematic resampling was selected was due to its low computational complexity in relation to the number of particles [17]. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 Multinomial Stratified Systematic
Figure 2.9:Three different uniform samples uk, generated from three
resam-pling techniques; multinomial resamresam-pling, stratified resamresam-pling and sys-tematic resampling.
2.2.5
UAV flight controller
The UAV controller addresses the task of using the estimated data from the filter to form a motion control for the seeker UAV.
One simple approach is to turn the seeker drone to the estimated bearing of the target drone and to stop if the estimated signal strength gets too high, to avoid a collision. A high signal strength implies that the seeker and target drone are close to each other.
Another method for calculating the motion control for a UAV is described in Sec-tion 2.2.6 where the Markov Decision Process is explained.
2.2.6
Optimality criterion
In this section different methods are described for finding optimality.
Monte Carlo simulations
Monte Carlo simulations consist of performing a repetitive simulation for a series of realisations. One simulation includes M realisations, all of them are indepen-dent and of length N . More accurately, it can be concentrated to the following formula ˆ A = 1 M M X i=1 xi[n], (2.30)
where x[n] = A + w[n] is a vector of length N and with mean A. w[n] is a white Gaussian noise of length N . This method requires that a new x[n] is computed for each realisation [19]. Monte Carlo simulations are useful in several aspects, in the formula above it is used to compute an estimator. Moreover, it can be used to evaluate the performance of a method [33]. However, the drawback is as stated previously, the method requires a new and independent realisation for each instance i. In some cases this can lead to expensive computations which can be avoided if the problem can be downsampled sufficiently [22]. In this project Monte Carlo simulations were used to evaluate the algorithm and to decide dif-ferent constants used in the algorithm.
Markov decision process (MDP)
A Markov Decisions Process, MDP, is a process of finding an optimal policy π∗
that minimises the cost function c(sk+∆k) for a time horizon T . The cost function
is a function of the state, where the new state sk+∆k is generated from the model
formulation using input uk. Thus each state is mapped to an action. Here uk is
an action from an action space U . The objective is to minimise the expected total cost, π∗(sk), to time horizon T such that [23]
π∗(sk) = argmin uk E T X t=0 c(sk+t∆k) . (2.31)
Here E denotes the expected value of the summed cost function c. Now follows an example from [22] where a cost function was designed to penalise the seeker drone if it got too close to the target. The implementation formulated a cost function as
c(sk) = λE n
||seeker pos. - target pos.|| < do (2.32)
where d was a pre-determined threshold. The condition within the expectation equalled 1 if its argument was true, otherwise 0. The weight λ decided the trade
off between the tracking and the collision avoidance. The higher λ was, the larger penalisation for near collisions. In this project the MDP was used as one of two different UAV controllers. A thoroughly description of how this automatic con-troller was implemented can be found in Section 3.4.1.
2.3
Noise distribution
In this section the noise distributions that were of interest for this project are presented.
2.3.1
Student-t distribution
The noise of the bearing is assumed to be Student-t distributed [20]. This distri-bution is quite similar to a Normal distridistri-bution, except that the probability of getting a value further away from the mean is higher. The probability density function is stated as follows
f (x|ν) = Γ( ν+1 2 ) Γ(ν2) 1 √ νπ 1 (1 +xν2)ν+12 (2.33) where v is a design parameter and denotes degrees of freedom. The gamma function, Γ , is as follows [7] Γ(n) = (n − 1)! = ∞ Z 0 tn−1etdt (2.34)
where n is any positive integer [6]. Figure 2.10, which is presented below, shows the probability density function of the Student-t distribution with two different degrees of freedom and a standard Normal distribution. Parameter v denotes the degrees of freedom for both the figure and the equation.
The student-t distribution was used in this project to simulate bearing noise. The noise was added to the calculated bearing, that originated from the target and the seeker drone’s exact positions (neglecting any GPS-error).
Section 3.5.2 contains a detailed explanation of the test framework setup for the bearing.
-5 0 5 Observation 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Probability Density
Student's t and Standard Normal pdfs Student's t Distr. with =3 Student's t Distr. with =1 Standard Normal Distr.
Figure 2.10:Two probability density functions of the Student-t distribution
and one probability density function of a Normal distribution.
2.3.2
Rice distribution
The noise of the signal strength is assumed to be Rice distributed, since the setup between the two drones assumed line of sight. This means that the signal has no or very low multi path propagation [32]. The probability density function for a Rice distribution is defined as follows
f (x|v, σ ) = x σ2e −(x2+v2) 2σ 2 I0 xv σ2 (2.35) where v is set according to
v = 1
r2. (2.36)
Here r is the distance between the drones and I0(z) is the modified Bessel
function of zeroth order [9]. The Bessel function In(z) is defined according to the
following equation In(z) = z n 2nn! " 1 − z 2 2(2n + 2)+ z4 2 · 4(2n + 2)(2n + 4)−... # (2.37) where n denotes the order of the function [5].
The probability density function of the Rice distribution is shown in Figure 2.11a and 2.11b. Figure 2.11a contains three probability density functions for differ-ent v but where σ is constant. Figure 2.11b contains three probability density functions for different σ but where v is constant.
0 2 4 6 8 Observation 0 0.1 0.2 0.3 0.4 0.5 0.6 Probability Density Rice pdf with = 1
Rice Distr. with =0 Rice Distr. with =2 Rice Distr. with =4
(a)σ = 1 where v is set to 0, 2 and 4.
0 2 4 6 8 Observation 0 0.1 0.2 0.3 0.4 0.5 Probability Density Rice pdf with = 1
Rice Distr. with =1 Rice Distr. with =2 Rice Distr. with =3
(b)v = 1 where σ is set to 1, 2 and 3.
Figure 2.11:Rice probability density function
The Rice distribution was used to simulate signal strength noise. The noise was added to the simulated signal strength, that originated from the target and the seeker drone’s exact positions (neglecting any GPS error).
Section 3.5.1 contains a detailed explanation of the test framework setup for the signal strength.
3
Method
In this chapter the methods used in order to answer the research questions are presented.
3.1
System overview
The project was divided in two parts, firstly all the code was tested and simulated on the computer, a theoretical test. The next step was to try the code in practice with drones. Each of these two parts was divided in three sub parts, namely simulation, tracking and heading adjustment. These three steps are explained in this chapter.
1. Theoretical test on computer • Simulation of signal • Tracking of position
• Heading adjustment of drone 2. Practical test with drones
• Simulation of signal • Tracking of position
• Heading adjustment of drone
In Section 1.4 it is mentioned that one of the limitations was that no antennas were used, neither for the theoretical or practical tests. Therefore the first step in the algorithm was to simulate the bearing and signal strength. When doing this theoretically the signal strength and bearing were calculated from the positions of the drones. To make the test more realistic, a significant amount of noise was
applied. In practice the drones built-in GPS position were used to get the position of the drones, thereafter the process was similar to the theoretical test, calculate the signal strength and bearing and apply a significant amount of noise. In this project it was assumed that the offset of the GPS systems could be neglected. The simulated noisy signal strength and bearing were sent to a tracking algo-rithm. The tracking algorithm was supposed to estimate the correct signal strength and bearing, in other words, to remove the noise. Two different tracking al-gorithms were tested and evaluated, one used a Particle filter and one used a Kalman filter. The implementation of the tracking algorithm is further discussed in Section 3.3.
The last step of the whole algorithm was the steering control. The task of the steering control was to change the direction and move the seeker drone to the direction of the target drone and to make sure the seeker drone took a smooth path. For each time the steering control algorithm was executed it generated a position the seeker drone was supposed to go to. In the theoretical test it gener-ated a position close to the seeker drone and in a perfect computer environment the simulated drone moved to that point exactly. The generated position was close to the seeker drone to mimic an actual movement. If the position would have been far away then it would not have been realistic, since the drone could not have reached the specific position in such a short amount of time. In practice, with real drones, it was not possible to send a position close to the seeker drone. The drone would think it already was at that position and stayed still. This was due to a margin in the drone’s GPS where it omitted positions which were within a radius of approximately five metres. Instead a position further away but in the right direction was sent, and before the drone reached that position a new posi-tion was sent, creating a smooth path.
For both the tracking algorithms two different steering controls were tested and evaluated. One steering control built on a Markov Decision Process and the other one, a more simple one, just turned to the estimated bearing.
A summary of the system overview is shown in Figure 3.1.
Simulation
Tracking
Algorithm
Steering
Control
Bearing and Signal Strength Particle Filter Kalman Filter Markov Decision process Turn to estimated bearing3.2
Setup
Two different setups for the estimation of the seeker and target positions were implemented. One of the setups used a global coordinate system only, whereas the other used both a global and a local coordinate system. In the end only the last setup was used.
3.2.1
Global coordinate system setup
This method made use of only one global coordinate system where the seeker and target positions were defined according to this system, see Figure 3.2.
α[n] α[n-1] Target at time step [n-1] Target at time step [n] Movement of seeker at time step [n-1] Global x axis Global y axis ϕ[n]
Figure 3.2:Illustration of the setup of the global coordinate system method.
In this setup two different angles were utilised in order for the seeker to estimate the position of the target, the relative bearing α and the world bearing φ. The world bearing φ was the angle between the seeker’s world east, or the x axis, and the direction to the target. The relative bearing was the angle between the pre-vious direction to the target and the current direction to the target, with respect to the seeker. In the first time step the relative bearing was assumed to be zero with respect to the world east, or x axis. The world bearing was calculated as a summation of all previously estimated relative bearings. The world bearing combined with the estimated signal strength were used to calculate the target position. The estimated relative bearing was used as an input to the automatic controllers. The cumulative property of the setup showed to be a weakness where
the error increased and when the standard deviation had a value of 20 the system completely broke down. More of this is found in Section 4.1.
3.2.2
Global and local coordinate system setup
This setup used two coordinate systems, one global and one local and can be seen in Figure 3.3. The global coordinate system did not change while the drones moved. The local coordinate system was normalised in the meaning that local y axis always followed the seeker drone’s direction. The relative bearing β was the angle between the local y axis and the target drone.
Global x axis Global y axis
Local x axis Local y axis
β
Figure 3.3: Illustration of the setup with a local and a global coordinate
system.
3.2.3
Conversion of GPS data to Cartesian coordinates
The positional data yielded from the GPS was the latitude, longitude and alti-tude, which can be summarised as World Geodetic System 1984 (WGS84) data. To obtain global coordinates this data had to be converted to a Cartesian coordi-nate system. The positional data was converted to Rikets koordinatsystem 1990 (RT90) data which is the coordinate system used for governmental maps in Swe-den. The coordinates were calculated according to the Gauss-Krüger projection with local mean meridians [24].
3.3
Tracking
For the tracking part different predefined paths for the target were used in order to check and verify the filters. As the seeker and/or target moved, the simulated
relative bearing and signal strength were calculated according to the new posi-tions. To each simulated observation an observation noise was added to mimic an authentic scenario. The simulated observation also yielded true data, which was the simulated data without noise. The true data was used for the evaluation of the filters.
3.3.1
Kalman Filter
The Kalman filter was implemented according to the summary in Chapter 2.2.3, the outline of the code is shown in the pseudo code below
Algorithm 2Scalar Kalman Filter for estimation of bearing
InputBearing signal β, scalar for amplification of filter m, variance σw,
vari-ance σu, scalar a
OutputEstimated bearing ˆβ
1: Initialise ˆβ 2: Initialise M 3: whileIncoming β do 4: β = a · ˆˆ β 5: M = √ a · M + (m · σw) 6: K = M/(σu+ M) 7: Update ˆβ → ˆβ + K · (β − ˆβ) 8: Update M → (1 − K)M 9: end while
The scalar m affected the Kalman filter in the way that a higher value would make the filter more prone to follow the noisy signal. A low value could, on the other hand, take time to stabilise. If the signal was not noisy or made quick changes it was a good idea to have a higher value on m, a noisier signal benefited from a lower value. To optimise the result of this project, Monte Carlo simulations were used to find the best value for m during different circumstances.
If there was any pre-known information about the situation the scalar a was used for this. In this project, when a seeker drone followed a target drone that took an unknown path, a was set to 1. There was no pre-known information.
For this project Algorithm 2 estimated the bearing β and the same code was used for the estimation of the signal strength. In the signal strength case the β was replaced with S and ˆβ with ˆS.
3.3.2
Particle Filter
The Particle filter was implemented according to the algorithm of the basic out-line of a Particle filter in Section 2.2.3.
As a first stage the noises were defined for the bearing and the signal strength, where their respective distributions were set according to the Student-T and the Rice distributions. The particles were initiated with a normalised weight where the sum of all weights equalled 1. The initial state for the bearing of each particle was Student-T distributed with 5 degrees of freedom. For the signal strength the initial state was Rice distributed where the design parameter v was set according to the estimated distance. When the distance was large v was set according to a large noise, when the distance was short v was set according to a small noise. The outline of the resampling made use of a cumulative summation for all of the particle weights at the current time step k, a cumulative distribution function (CDF). A temporary vector representing the weights of the particles was used and compared to the CDF. The temporary vector consisted of values in increasing or-der from 0 to 1, where the last position held a weight with a value close to 1. The outline of the systematic resampling method is described below. The resampled
Algorithm 3Systematic resampling
InputParticle weights wk, model state sk
OutputResampled model state sjres, resampled particle weights wresj
1: Create CDF c
2: Create temporary vector u
3: Set u0∼U [0,N1s]
4: Set i = 0
5: for all particles j do
6: uj= u0+ (j − 1)N1s 7: whileuj> ci do 8: i = i + 1 9: end while 10: Assign sjres= si 11: Assign wresj =N1 s 12: end for
states sresk were the states which were used in the Particle filter for future time
steps k + ∆. The relation between the CDF and the particle weights is illustrated in Figure 3.4.
3.4
Flight control
For the flight controlling two different methods were implemented. One of them simply moved according to the relative bearing. The other one was based on a Markov Decision Process where different control inputs were explored. The inputs to the controllers were the estimated relative bearing and signal strength obtained from the tracking filters.
w0 w1 ... wN ... CDF Weight 0 1 Probability Particle weights
Figure 3.4:Generation of a cumulative sum using the weights from a particle
3.4.1
MDP controller
For the MDP controller eight different inputs were explored for a time horizon of four time steps. The resulting optimal policy from the MDP selected the control input which yielded the smallest cost, given a certain state. The cost function used for this purpose was stated as follows
c(sk) = ||target pos. in next time step - seeker pos||. (3.1)
The optimal policy was designed according to the following equation
π∗(sk) = argmin uk T X t=0 c(sk+t∆k). (3.2)
The seeker position was updated in a loop for each new propagation towards the time horizon, given one of the eight control inputs. The control inputs consisted
of varying the bearing in the interval [−20◦
, 20◦
] with respect to the estimated relative bearing at time step k. The target position in Equation 3.1 was the pre-dicted target position for the next time step k + 1 given the last and second to last estimated target positions. A linear equation for these two previous positions was calculated in order to obtain the predicted position. The seeker position in the same equation was the current seeker position for the current time index t. The seeker position was then updated after the optimal policy and the collision handling.
3.4.2
Relative bearing controller
Compared to the MDP controller, this method was much simpler. With this method the seeker moved according to the estimated relative bearing. In prac-tice this meant that the seeker drone rotated and headed toward the target drone according to the latest estimated relative bearing.
3.4.3
Collision handling
The collision detection and handling commenced right after the automatic con-trolling process ended. Firstly the detection checked if the distance between the seeker and the target, given the chosen steering instruction from the automatic
controlling process, exceeded threshold d1and d2. If that was the case then the
seeker should act according to the steering decision. If the distance between
seeker and target exceeded d1, but did not exceed d2, then the chosen steering
should not be acted upon. Instead the seeker should then remain at the same position as before. If the chosen steering instruction resulted in a distance which
did not exceed threshold d2nor d1, then the seeker should move backwards. An
Target
d1
d2
Figure 3.5:The thresholds d1and d2for the collision handling of seeker and
target.
3.5
Test framework setup
In both the theoretical and practical tests the bearing and signal strength were estimated. For both cases the systems were tested in different settings. What hap-pens if the outgoing signal from the target drone is weak? What haphap-pens if the (imaginary) antennas on the seeker drone work poorly? Together with experts at the Swedish Defence Research Agency (FOI) we generated two models to simu-late the problems above. The model for the signal strength is discussed in Section 3.5.1 and the model for the bearing is discussed in Section 3.5.2.
3.5.1
Test model for signal strength
As mentioned in Chapter 2.3.2 the noise distribution for the signal strength was assumed to be Rice distributed. Different amounts of noise were used to simulate different amount of outgoing signal strength from the target drone. The amount of noise on the signal also depended on the distance between the drones. For both the problems Rician fading was used. The Rician fading used the parameter K to illustrate the ratio between the power of the LOS component and the other scattered paths. Different intervals of K were used to simulate different amount of outgoing signal strength from the target drone. The relation between K and
the variance of the noise σ2, is shown in the equations below. The variance was
calculated according to the following equation [29] σ2= v
2
2K. (3.3)
v = 1
r2. (3.4)
Here r was the distance between the two drones. The K values could then be calculated accordingly,
K = v
2
2σ2. (3.5)
The variance was then used in the Rice distribution to calculate the noise for the signal amplitude as shown in Equation 2.35.
In this project the parameter K depended linearly on the distance between the drones. A smaller distance generated a higher value on K and therefore a smaller amount of noise on the signal strength. The closest distance possible, namely 0 meter, generated a K value that was 10 times higher compared to the greatest distance. For this project the limit was 100 meters because of the wifi-range. When discussing this with our advisers at FOI we decided to test K-intervals from 1 to 10 for strong noise and from 100 to 1000 for small noise.
Table 3.1: Summation of different K-intervals and how they relate to the
signal strength noise.
K-interval
Small noise 100 - 1000
Medium noise 10 - 100
Large noise 1 - 10
Table 3.1 shows a summation of how different K intervals were used to simu-late different signal strength noises. Figure 3.6 shows how the K value changed linearly depending on the distance between the two drones. The K value is shown for three different K-intervals.
To verify that the simulated noise behaved as expected, the error in distance was examined in two different ways. The error in distance was the subtraction between two signals, namely the signal with applied noise and the same signal before applying the noise, the ground truth. This gave a geometric error in me-ters of how far the estimated position was from the exact position. Both tests were performed on a computer and the two drones were simulated.
During the first test the seeker drone followed the target drone three times. In all three times the target drone flew the same path and the seeker started from the same position. It took the target drone a little over 350 samples to finish its path. For this test setup there was no K interval, the value of K was constant and did not depend on the distance between the drones. Three different values on K were used, K = 1, K = 100 and K = 1000. Figure 3.7 shows the error in distance for these three tests. From this it is possible to see that the error decreases when the K increases.
0 20 40 60 80 100 Distance r [m] 0 200 400 600 800 1000 K-va lue
K-value depening on distance
K-interval 1-10 K-interval 10-100 K-interval 100-1000
Figure 3.6:Shows how the K-value depended linearly on the distance r
be-tween the two drones for three different noise settings.
0 50 100 150 200 250 300 350 Sample 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Error in distance
Error in distance for different values K
K = 1 K = 100 K = 1000
[m]
Figure 3.7:Error in distance between signal with noise and the signal
with-out noise and three different values on K.
The goal of the second simulation test was to verify that the error decreased when the distance between the drone decreased. For this the K value changed linearly depending on distance between 10 and 100. This K interval is shown in Figure 3.6 as the orange line. The simulated drones started approximately 50 meters from each other, after 25 samples the distance was approximately 10 meters and after hardly the 50 samples the distance was approximately 15 meters. Figure
36 3 Method 0 10 20 30 40 50 0.5 1.0 1.5 2.0 2.5 Err or in dis tan ce [m ] 0 10 20 30 40 50 Sample 10 20 30 40 50 Dis tan ce [m ]
Distance between the dronesSample
0 10 20 30 40 50 0.5 1.0 1.5 2.0 2.5 3.0 Err or in dis tan ce [m ] Error in distance 0 10 20 30 40 50 Sample 10 20 30 40 50 Dis tan ce [m ]
Distance between the drones
(a) 0 10 20 30 40 50 0.5 1.0 1.5 2.0 2.5 3.0 Err or in dis tan ce [m ] Error in distance 0 10 20 30 40 50 Sample 10 20 30 40 50 Dis tan ce [m ]
Distance between the dronesSample
0 10 20 30 40 50 0.5 1.0 1.5 2.0 2.5 3.0 Err or in dis tan ce [m ] Error in distance 0 10 20 30 40 50 Sample 10 20 30 40 50 Dis tan ce [m ]
Distance between the drones (b)
Figure 3.8: Two plots of a simulated test of ∼ 50 samples where the seeker
drone approached the target drone.
(a): Plot of how the distance between the drones changed. (b): Plot of how the error in distance changed.
3.8(a) shows the distance between the drones for the described simulation. Fig-ure 3.8(b) shows the error in distance generated from the simulation. Likewise as before, the error in distance is the subtraction between two signals, namely the signal with applied noise and the same signal before applying the noise, the ground truth. In the figure it is possible to see that the error in distance reduced from ∼ 3 meter to ∼ 0.5 meter depending on the distance of the drones. Impor-tant to mention is that the error not only depended on the value K but also on the signal strength, that grew stronger when the drones moved closer.
From these two examinations the conclusion was drawn that the error depended on different K intervals (Figure 3.7) and that the error also depended on the dis-tance between the drones (Figure 3.8).
