Adaptive Filtration of Parameters of the UAV Movement Based on the TDOA-Measurement Sensor Networks

ABSTRACT:

In modern conditions unmanned aerial vehicles (UAVs) generate new classes of threats, including their use for terrorist purposes. A feature of modern UAVs is the ability to perform sudden maneuvers and to keep the same position in the point in space. For the description of the UAV movement with various types of maneuver it is used a rectangular coordinate system. We use the model in the form of stochastic dynamic system with random structure in the discrete time in which the change type UAV movement occurs at random times. When a UAV emits a sign, its location can be determined by wireless sensor networks (WSN) using the TDOA method. On the basis of a mathematical apparatus of the mixed Markov processes for in discrete time optimal and quasi-optimal adaptive algorithms for filtration of UAV movement parameters based on the TDOA-measurement, sensor networks are synthesized. Devices that realize these algorithms are multichannel and belong to the class of devices with feedback between channels. At the same time, in a quasi-optimal algorithm, a sequential procedure of the arriving measurements from sensors of a sensor network is realized, which allows to avoid the inversion of large-dimensional matrices. An analysis of the quasi-optimal adaptive algorithm is performed using statistical modeling. On the intervals of hovering and of the UAV movements without maneuver, the developed algorithm allows to increase significantly the accuracy of the estimation of the UAV coordinates, and also to recognize various types of its movement with high probability level.

KEYWORDS:
UAV; Wireless sensor networks; Optimal and quasi-optimal adaptive algorithms; Parameters of the movement; TDOA

INTRODUCTION

Unmanned aerial vehicles (UAVs) are beginning to play an increasingly important role. The main arguments for their dissemination are low cost, high mobility and the ability to transfer photo-video information in real time. On the other hand, the availability and massive use of UAVs has led to a new class of threats (Nonami et al. 2010; Wallace and Loffi 2015; Sathyamoorthy 2015; Solodov et al. 2018; Friese et al. 2016; Card 2018; Dronebouncer 2018).

The development of wireless technology and microelectronics has led to the emergence and widespread use of wireless sensor networks (WSN), which are distributed networks self-organizing from miniature sensors that exchange information through wireless technology. One of the main uses of WSN is the creation of various monitoring and tracking systems of moving radio sources, including UAVs (Gemayel et al. 2014; Wan et al. 2016).

To determine the UAV location using wireless sensor networks, TDOA measurements are widely used (Rullan-Lara et al. 2013; Makki et al. 2016). At TDOA measurement, synchronization between WSN and a radio source is not necessary, which allows its use for positioning the unknown UAV.

The UAV trajectory as a rule consists of non-overlapping sections with different types of movement. It can be presented in the form of stochastic process from which probability characteristics change at random instants of time. A convenient mathematical model of such processes are stochastic discrete dynamic systems with a random structure, which are adequate to solvable tasks when realization of algorithms on digital computers. To increase the accuracy of the UAV coordinates estimate on sections with different types of movement, various measurement processing algorithms must be used. However, the type of movement is, as a rule, unknown. Therefore, besides estimating the UAV coordinates, it is also necessary to solve the task of recognizing its motion type.

In the work of Tovkach and Zhuk (2017a) the adaptive algorithm is synthesized, with estimation of the UAV movement parameters and its movement type. Observations at synthesis of the algorithm are the UAV location coordinates in the rectangular coordinate system, which are determined on the basis of TDOA measurements using the least squares method. The disadvantage of the obtained algorithm is that definition the UAV location is carried out after arrival of measurements from all sensors, which leads to the necessity inversion of large-dimensional matrices.

Therefore, the synthesis of adaptive algorithms for the UAV movement parameters estimate based on the TDOA-measurement sensor networks, which provide serial processing of measurements, is important and do not require the inversion of large-dimensional matrices.

PROBLEM DEFINITION

The UAV movement with various types of maneuver in rectangular system of coordinates can be described by a stochastic dynamic system with random structure (Eq. 1) (Zhuk 1989):

(1) ${\mathit{u}}^{\mathit{tr}}\left(\mathit{k}\right)={\mathit{F}}_{\mathit{j}}^{\mathit{tr}}{\mathit{u}}^{\mathit{tr}}\left(\mathit{k}-1\right)+{\mathit{G}}_{\mathit{j}}^{\mathit{tr}}{\mathit{\omega }}^{\mathit{tr}}\left(\mathit{k}\right),\mathit{j}=\stackrel{─}{1,\mathit{M}},$

where u^tr (k) is the state vector; F_j^tr and G_j^tr are the matrixes describing various types of movement; ω^tr(k) is the uncorrelated sequence of Gaussian vectors with a single correlation matrix.

During propagation of a signal from the UAV to sensors of network, its coordinates do not change. Therefore, it is assumed that the measurements of the differences of distances between the sensors and the reference sensor are received in the k-th moment of time. The moment of receipt of a signal to the basic sensor with coordinates of x⁰, y⁰, z⁰ is used as the k-th moment of time.

Distance difference measurement errors are correlated because they contain reference sensor measurement error. The presence of correlated errors makes it difficult to use traditional recurrent target coordinates evaluation algorithms. This difficulty can be avoided through the introduction of error υ⁰(k) into the state vector of the estimated parameters (Tovkach and Zhuk 2017b), which can be represented in the form u^T(k) = (u^trT(k), υ⁰(k)). In this case, the model of the process to be filtered is transformed to Eq. 2:

(2) $\mathit{u}\left(\mathit{k}\right)={\mathit{F}}_{\mathit{j}}\mathit{u}\left(\mathit{k}-1\right)+{\mathit{G}}_{\mathit{j}}\mathit{\omega }\left(\mathit{k}\right),\mathit{j}=\stackrel{─}{1,\mathit{M}},$

where ω^T(k) = (ω^trT(k), ω_Δ(k)); ω_Δ(k) is the uncorrelated sequence of Gaussian values with a single dispersion;

${\mathit{F}}_{\mathit{j}}=\left[\begin{array}{cc}{\mathit{F}}_{\mathit{j}}^{\mathit{tr}}& 0\\ 0& 0\end{array}\right];{\mathit{G}}_{\mathit{j}}=\left[\begin{array}{cc}{\mathit{G}}_{\mathit{j}}^{\mathit{tr}}& 0\\ 0& {\mathit{\sigma }}_{∆}\end{array}\right].$

In this case, the observation equation that describes the process of measuring the coordinates of UAV by sensor network is Eq. 3:

(3) ${∆}^{\mathit{s}}\left(\mathit{k}\right)={\mathit{h}}^{\mathit{s}}\left(\mathit{u}\left(\mathit{k}\right)\right)+{\mathit{\upsilon }}^{\mathit{s}}\left(\mathit{k}\right),\mathit{s}=\stackrel{─}{1,\mathit{S}},$

where Δ^s(k) is the measured difference of distances between the s-th sensor and the reference sensor in the k-th moment of time: s = 1,S; υ^s (k) is the measurement error of s-th sensor with dispersion d^s; h^s(u) is the nonlinear function, described by Eq. 4:

(4) ${\mathit{h}}^{\mathit{s}}\left(\mathit{u}\left(\mathit{k}\right)\right)=\sqrt{{\left(\mathit{x}\left(\mathit{k}\right)-{\mathit{x}}_{\mathit{s}}\right)}^2+{\left(\mathit{y}\left(\mathit{k}\right)-{\mathit{y}}_{\mathit{s}}\right)}^2+{\left(\mathit{z}\left(\mathit{k}\right)-{\mathit{z}}_{\mathit{s}}\right)}^2}-\sqrt{{\mathit{x}}^2\left(\mathit{k}\right)+{\mathit{y}}^2\left(\mathit{k}\right)+{\mathit{z}}^2\left(\mathit{k}\right)}-{\mathit{\upsilon }}^0\left(\mathit{k}\right);$

where x(k), y(k), z(k) are coordinates of the UAV position; x^s, y^s, z^s are coordinates of the s-th sensor position; υ⁰(k) is the error of the distances measurement by the reference sensor with dispersion d⁰. Measurements obtained from all sensors are denoted as a vector u_Δ(k) = (Δ¹(k),…,Δ^S(k))^T.

On the basis of Eqs. 2 and 3 it is necessary to synthesize optimal and quasi-optimal adaptive filtration algorithms of UAV movement parameters according to sensor network, in which the detection of various types of motion is carried out.

DEVELOPMENT

Equation 2 describes the process to be estimated, and Eq. 3 the formation process of data accessible to observation. For a description of the type of maneuver there is used the random variable a^j(k), j = 1,M which belongs to chains of Markov, accepting the M values with a matrix of probabilities of transitions Π_ij (k, k - 1) i, j = 1,M and initial probabilities p_i (0), i = 1,M. Measurement errors υ⁵ (k), s = 1,S are independent.

Applying the technique of a condition vector expansion (Zhuk 1988) it is possible to show that the expanded process (u(k), a_j(k)) belongs to the class of the mixed Markov processes in discrete time with transition probability density (PD) (Eq. 5):

(5) $\Pi \left(\mathit{u}\left(\mathit{k}\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}-1\right),{\mathit{a}}_{\mathit{i}}\left(\mathit{k}\right)\right)=\Pi \left(\mathit{u}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}-1\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right){\Pi }_{\mathit{ij}}\left(\mathit{k},\mathit{k}-1\right),$

where Π(u(k) | u(k - 1), a_j(k)) is the conditional PD defined by Eq. 2.

The problem of synthesis of an optimal algorithm comes down to calculation of aposteriori PD of expanded process W(u(k), a_j(k)) = P(u(k), a_j(k) | U_Δ(k)), where U _Δ(k) = u _Δ(1),…,u _Δ(k) are the sequences of measurements received to till k-th moment, included. Introducing also designation of the extrapolated PD P(u(k), a_j (k) | U_Δ(k - 1)) for the expanded process and following the technique of synthesis (Zhuk 1989), the optimal algorithm of a filtration can be presented in the form of two the recurrent Eqs. 6 and 7:

(6) ${\mathit{W}}^{*}\left(\mathit{u}\left(\mathit{k}\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right)=\sum _{\mathit{i}=1}^{\mathit{M}}{\Pi }_{\mathit{ij}}\left(\mathit{k},\mathit{k}-1\right)\times \underset{-\infty }{\overset{\infty }{\int }}\Pi \left(\mathit{u}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}-1\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right)\times \mathit{W}\left(\mathit{u}\left(\mathit{k}-1\right),{\mathit{a}}_{\mathit{i}}\left(\mathit{k}-1\right)\right)\mathit{du}\left(\mathit{k}-1\right);$

(7) $\mathit{W}\left(\mathit{u}\left(\mathit{k}\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right)=\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}\right)\right)\times {\mathit{W}}^{*}\left(\mathit{u}\left(\mathit{k}\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right)/\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right),$

where $\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}\right)\right)=\prod _{\mathit{s}=1}^{\mathit{s}}\mathit{P}\left({∆}^{\mathit{s}}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}\right)\right)$ is one-step-likelihood function, defined on the basis of measurements Δ ^s(k), s = 1,S; P(Δ^s(k) | u(k)) is the conditional PD, defined by Eq. 6; P(u_Δ(k) | U_Δ(k - 1)) is the conditional PD determined by:

$\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right)=\sum _{\mathit{j}=1}^{\mathit{M}}\underset{-\infty }{\overset{\infty }{\int }}\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}\right)\right){\mathit{W}}^{*}\left(\mathit{u}\left(\mathit{k}\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right)\mathit{du}\left(\mathit{k}\right).$

Equation 6 describes evolution of the extrapolated PD W^*(u(k), a_j(k)) and is also an optimal algorithm of extrapolation of mixed Markov processes (u(k), a_j(k)) for one step. With the help of Eq. 7 the clarifying of extrapolated PD based on the obtained measurement u_Δ(k) is made and a posteriori PD W(u(k), a_j(k)) is defined.

Further transformation of Eqs. 6 and 7 can be executed by means of the theorem of multiplication of probabilities. In this case the optimal filtering algorithm can be represented as the following system of recurrent Eqs. 8 to 11:

(8) ${\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{k}\right)=\sum _{\mathit{i}=1}^{\mathit{M}}{\Pi }_{\mathit{ij}}\left(\mathit{k},\mathit{k}-1\right){\mathit{W}}_{\mathit{i}}\left(\mathit{k}-1\right);$

(9) ${\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{u}\left(\mathit{k}\right)\right)=\sum _{\mathit{i}=1}^{\mathit{M}}{\Pi }_{\mathit{ij}}\left(\mathit{k},\mathit{k}-1\right){\mathit{W}}_{\mathit{i}}\left(\mathit{k}-1\right)\times \underset{-\infty }{\overset{\infty }{\int }}\Pi \left(\mathit{u}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}-1\right),{\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right)\right)\times {\mathit{W}}_{\mathit{i}}\left(\mathit{u}\left(\mathit{k}-1\right)\right)\mathit{du}\left(\mathit{k}-1\right)/{\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{k}\right);$

(10) ${\mathit{W}}_{\mathit{j}}\left(\mathit{u}\left(\mathit{k}\right)\right)=\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}\right)\right)\times {\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{u}\left(\mathit{k}\right)\right)/\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right),{\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right);$

(11) ${\mathit{W}}_{\mathit{j}}\left(\mathit{k}\right)=\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right),{\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right)\times {\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{k}\right)/\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right),$

where W^*_j(u(k)) = P(u(k) | a_j(k), U_Δ(k - 1)) and W_j(u(k)) = P(u(k) | a_j(k), U_Δ(k)) are the conditional extrapolated and a posteriori PD of the vector u(k) on condition a_j(k); W^*_j(k) = P(a_j(k), U_Δ(k - 1)) and W_j(k) = P(a_j(k), U_Δ(k)) are extrapolated and a posteriori probabilities a_j(k); P(u(k) | a_j(k), U_Δ(k - 1)) is the conditional PD determined by:

$\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right),{\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right)=\underset{-\infty }{\overset{\infty }{\int }}\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid \mathit{u}\left(\mathit{k}\right)\right){\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{u}\left(\mathit{k}\right)\right)\mathit{du}\left(\mathit{k}\right),\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right),$

determined by: $\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right)=\sum _{\mathit{j}=1}^{\mathit{M}}\mathit{P}\left({\mathit{u}}_{∆}\left(\mathit{k}\right)\mid {\mathit{a}}_{\mathit{j}}\left(\mathit{k}\right),{\mathit{U}}_{∆}\left(\mathit{k}-1\right)\right){\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{k}\right)$

The initial conditions for the algorithm Eqs. 8 to 11 have the form W_i(0) = p_i(0), W_i(u(0)) = P(u(0)), i = 1,M.

Using the Eqs. 8 and 11, extrapolated W^*_j(k) and a posteriori W_j(k) probabilities are calculated and Eqs. 9 and 10 are conditional extrapolated W^*_j(u(k)) and a posteriori W_j(k) PD. The optimal adaptive filter (Eqs. 8 to 11) contains M channels and belongs to the class of devices with feedbacks between channels (Zhuk 1988; 1989).

Large computational costs are required for accomplishing the optimal algorithm. When synthesizing a quasi-optimal algorithm, there will be used a sequential method of processing incoming data. There will enter the index l, which characterizes the sequence of receipts of measurements from the corresponding sensors of network. The quasi-optimal algorithm of an adaptive filtration can be received by decomposition of vector-valued function (Eq. 3), in a Taylor series in the vicinity of the points u(k) and restriction by linear members of decomposition, and also application on each step of Gaussian approximation of the conditional extrapolated PD W^*_j(u(k)). In this case the equation of calculation of the conditional extrapolated PD W*_j(u(k)) (Eq. 9) comes down to calculation of its first u^*_j(k) and second P^*_j(k) moments by Eqs. 12 and 13 (Zhuk 1989):

(12) ${\mathit{u}}_{\mathit{j}}^{*}\left(\mathit{k}\right)=\sum _{\mathit{i}=1}^{\mathit{M}}{\Pi }_{\mathit{ij}}\left(\mathit{k},\mathit{k}-1\right){\mathit{W}}_{\mathit{i}}\left(\mathit{k}-1\right)\mathit{F}{û}_{\mathit{i}}\left(\mathit{k}-1\right)/{\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{k}\right);$

(13) ${\mathit{P}}_{\mathit{j}}^{*}\left(\mathit{k}\right)=\sum _{\mathit{i}=1}^{\mathit{M}}{\Pi }_{\mathit{ij}}\left(\mathit{k},\mathit{k}-1\right){\mathit{W}}_{\mathit{i}}\left(\mathit{k}-1\right)\times \left\{{\mathit{F}}_{\mathit{j}}{\hat{\mathit{P}}}_{\mathit{i}}\left(\mathit{k}-1\right){\mathit{F}}_{\mathit{j}}^{\mathit{T}}+{\mathit{G}}_{\mathit{j}}{\mathit{G}}_{\mathit{j}}^{\mathit{T}}\right\}/{\mathit{W}}_{\mathit{j}}^{*}\left(\mathit{k}\right).$

Calculation of conditional a posteriori PD W_j(u(k)) (Eq. 10) at sequentially processing of the arriving measurements Δ^l(k), l = 1,S comes down to calculation of its first û_j(k) and second P̂_j(k) moments using the recurrent procedure (Eqs. 14 to 16) (Evlanov and Zhuk 1990):

(14) ${\mathit{K}}_{\mathit{j}}^{\mathit{l}}\left(\mathit{k}\right)={\hat{\mathit{P}}}_{\mathit{J}}^{\mathit{l}-1}\left(\mathit{k}\right)\frac{\partial {\mathit{h}}^{\mathit{lT}}\left({û}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)\right)}{\partial \mathit{u}\left(\mathit{k}\right)}\left({\mathit{D}}_{\mathit{j}}^{\mathit{l}}\left(\mathit{k}\right)\right){}^{-1};$

(15) ${û}_{\mathit{j}}^{\mathit{l}}={û}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)+{\mathit{K}}_{\mathit{j}}^{\mathit{l}}\left(\mathit{k}\right)\left({∆}^{\mathit{l}}\left(\mathit{k}\right)-{\mathit{h}}^{\mathit{l}}\left({û}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)\right)\right);$

(16) ${\hat{\mathit{P}}}_{\mathit{j}}^{\mathit{l}}\left(\mathit{k}\right)={\hat{\mathit{P}}}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{K}}_{\mathit{j}}^{\mathit{l}}\left(\mathit{k}\right)\frac{\partial {\mathit{h}}^{\mathit{l}}\left({û}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)\right)}{\partial \mathit{u}\left(\mathit{k}\right)}{\hat{\mathit{P}}}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right),$

where ûl_j(k), P̂l_j(k) are the expected value and correlation matrix of conditional a posteriori PD W_j(u(k)), refined by measurements Δ^s(k), s = 1, l. In this case the initial conditions for the procedure Eqs. 14 to 16 at l = 0 has the form û⁰_j = û^*_j(k), P̂⁰_j(k) = P̂^*_j(k), j = 1, M, and û_j(k) = û^L_j(k), P̂_j(k), j = 1, M. The advantage of this method is the possibility to organize processing with limited computing resources. It is especially convenient in case of nonsynchronous moments of input data arrival in comparison with algorithms of parallel processing (Evlanov and Zhuk 1990; Tovkach and Zhuk 2017c).

The algorithm of the filtration of discrete components does not change and is described by the Eqs. 8 and 11.

In this case conditional PD P(u_Δ(k) | a_j(k), U_Δ(k - 1)) = N(H(u^*_j(k)), D_j(k)) is Gaussian, and correlation matrix D_j(k) is determined by Eq. 17:

(17) ${\mathit{D}}_{\mathit{j}}^{\mathit{l}}\left(\mathit{k}\right)=\frac{\partial {\mathit{h}}^{\mathit{l}}\left({û}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)\right)}{\partial \mathit{u}\left(\mathit{k}\right)}{\hat{\mathit{P}}}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)\frac{\partial {\mathit{h}}^{\mathit{lT}}\left({û}_{\mathit{j}}^{\mathit{l}-1}\left(\mathit{k}\right)\right)}{\partial \mathit{u}\left(\mathit{k}\right)}+{\mathit{d}}^{\mathit{l}}\left(\mathit{k}\right).$

The quasi-optimal device that implements the algorithm Eqs. 8, 11 to 16 is also multichannel with number of channels M and generally keeps the structure of the optimal device. Their difference is that in quasi-optimal device only the first and second moments of conditional PD W^*_j(u(k)) and W_j(u(k)) are calculated. In this case the quasi-optimal algorithm allows keeping representation of a posteriori PD W_j(u(k)) in the form of the sum of the M Gaussian densities.

The block diagram of the j-th channel of such quasi-optimal filter is shown in Fig. 1.

Figure 1
Block diagram of one channel quasi-optimal filter at consecutive calculations.

THE EFFECTIVENESS OF THE ALGORITHM ANALYSIS

The operability of the obtained quasi-optimal adaptive algorithm of the UAV movement parameters estimate on the basis of TDOA measurements Eqs. 8, 11 to 16 is checked by the Monte Carlo method.

The sensor network (Fig. 2) consists of nine sensors with coordinates S₀(0; 0; 0), S₁(0; 100; 44), S₂(70.71; 70.71; 6), S₃(100; 0; 48), S₄(70.71; -70.71; 10), S₅(0; -100; 52), S₆(-70.71; -70.71; 14), S₇(-100; 0; 56), and S₈(-70.71; 70.71; 18).

Figure 2
Configuration of the sensor network with nine sensors and the trajectory of UAV movement.

The created test trajectory of the movement UAV (Fig. 1) consists of nine sections: first 1 ≤ k ≤ 19 - uniform motion; second 20 ≤ k ≤ 27 - maneuver, height reduction; third 28 ≤ k ≤ 42- uniform motion; fourth 43 ≤ k ≤ 45 - maneuver; fifth 46 ≤ k ≤ 59 - uniform motion; sixth 60 ≤ k ≤ 67 - maneuver, rise to height; seventh 68 ≤ k ≤ 85 - uniform motion; eighth 85 ≤ k ≤ 114 - hanging; ninth 115 ≤ k ≤ 130 - uniform motion. Error RMS of measurement sensors of sensor network σ_v = 2.4 m, step of sampling Т = 1 s. The simulation was carried out in a hundred realizations with various measurement errors.

A model with a random structure (Eq. 2) takes into account the three main types of UAV movement M = 3: hanging j = 1, almost uniform motion j = 2, movement with maneuver j = 3.

The state vector u^T(k) = (x(k), ẋ(k), ẍ(k), y(k), ẏ(k), ÿ(k), z(k), ż(k), z̈(k), v⁰(k)) includes position coordinates, velocity and acceleration along the axes X, Y, Z. The matrixes F_j(k, k - 1), G_j(k), j = 1,3 have the form:

${\mathit{F}}_{\mathit{j}}^{\mathit{tr}}\left(\mathit{k},\mathit{k}-1\right)=\left[\begin{array}{ccc}{\mathit{F}}_{\mathit{j}}^{\mathit{b}}& 0& 0\\ 0& {\mathit{F}}_{\mathit{j}}^{\mathit{b}}& 0\\ 0& 0& {\mathit{F}}_{\mathit{j}}^{\mathit{b}}\end{array}\right],{\mathit{G}}_{\mathit{j}}^{\mathit{tr}}\left(\mathit{k}\right)=\left[\begin{array}{ccc}{\mathit{G}}_{\mathit{j}}^{\mathit{b}}& 0& 0\\ 0& {\mathit{G}}_{\mathit{j}}^{\mathit{b}}& 0\\ 0& 0& {\mathit{G}}_{\mathit{j}}^{\mathit{b}}\end{array}\right],$

where F^b_j, G^b_j and j = 1,3 have the form

${\mathit{F}}_1^{\mathit{b}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right];{\mathit{F}}_2^{\mathit{b}}=\left[\begin{array}{ccc}1& \mathit{T}& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right];{\mathit{F}}_3^{\mathit{b}}=\left[\begin{array}{ccc}1& \mathit{T}& \frac{{\mathit{T}}^2}{2}\\ 0& 1& \mathit{T}\\ 0& 0& 1\end{array}\right];{\mathit{G}}_1^{\mathit{b}}=\left[\begin{array}{c}{\mathit{a}}_1·\mathit{T}\\ 0\\ 0\end{array}\right];{\mathit{G}}_2^{\mathit{b}}=\left[\begin{array}{c}\frac{{\mathit{a}}_2·{\mathit{T}}^2}{2}\\ {\mathit{a}}_2·\mathit{T}\\ 0\end{array}\right];{\mathit{G}}_3^{\mathit{b}}=\left[\begin{array}{c}\frac{{\mathit{a}}_3·{\mathit{T}}^3}{6}\\ \frac{{\mathit{a}}_3·{\mathit{T}}^2}{2}\\ {\mathit{a}}_3·\mathit{T}\end{array}\right];{\mathit{a}}_1,{\mathit{a}}_2,{\mathit{a}}_3$

are RMS of random fluctuations of speed, acceleration and speed of change of acceleration of the UAV, respectively. When modeling relied: a₁ = 0.05 m/s; a₂ = 0.01 m/s²; a₃ = 6 m/s³.

The partial derivatives ∂h^l(û_j^l-1(k))/∂u(k) included in Eqs. 14 to 17 have the form:

$\frac{\partial {\mathit{h}}^{\mathit{l}}\left({\hat{\mathit{u}}}_{\mathit{j}}^{-\mathit{l}1}\left(\mathit{k}\right)\right)}{\partial \mathit{u}\left(\mathit{k}\right)}=\left[\frac{{\hat{\mathit{x}}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{x}}_{\mathit{s}}}{{\hat{\mathit{R}}}^{\mathit{l}}}-\frac{{\hat{\mathit{x}}}^{\mathit{l}-1}\left(\mathit{k}\right)}{{\hat{\mathit{R}}}^0};0;0;\frac{{\hat{\mathit{y}}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{y}}_{\mathit{s}}}{{\hat{\mathit{R}}}^{\mathit{l}}}-\frac{{\hat{\mathit{y}}}^{\mathit{l}-2}\left(\mathit{k}\right)}{{\hat{\mathit{R}}}^0};0;0;\frac{{\hat{\mathit{z}}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{z}}_{\mathit{s}}}{{\hat{\mathit{R}}}^{\mathit{l}}}-\frac{{\hat{\mathit{z}}}^{\mathit{l}-2}\left(\mathit{k}\right)}{{\hat{\mathit{R}}}^0};0;0;-1\right],$

where ${\hat{\mathit{R}}}^{\mathit{l}}=\sqrt{{\left({\hat{\mathit{x}}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{x}}_{\mathit{s}}\right)}^2+{\left({\hat{\mathit{y}}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{y}}_{\mathit{s}}\right)}^2+{\left({\hat{\mathit{z}}}^{\mathit{l}-1}\left(\mathit{k}\right)-{\mathit{z}}_{\mathit{s}}\right)}^2}$

estimates the distance between l-th sensor and UAV, l = 1,S; ${\hat{\mathit{R}}}^0=\sqrt{{\hat{\mathit{x}}}^{\mathit{l}-1}{\left(\mathit{k}\right)}^2+{\hat{\mathit{y}}}^{\mathit{l}-1}{\left(\mathit{k}\right)}^2+{\hat{\mathit{z}}}^{\mathit{l}-1}{\left(\mathit{k}\right)}^2}$

estimates the distance between the reference sensor and the UAV.

The initial conditions û_j(0), P̂_j(0), j = 1,3 for the first hypothesis j = 1 were created on the current measurements, and for j = 2,3 according to the observations at the two neighboring steps.

Figure 3 shows the dependences of recognition probability of movement of the first, second and third types, respectively, curves 1, 2, 3, obtained by statistical modeling. The probability of recognizing the hovering and almost uniform motion of the UAV is close to unity by the algorithm obtained. The probability of recognizing the maneuver is lower due to its short duration.

Figure 3
The probability of determining maneuvers.

Figure 4 shows dependences of expected value (curve 1) and RMS (curve 2) errors of the position estimation of the UAV along coordinates of X, Y, Z and also RMS (curve 3) errors of assessment calculated by the adaptive algorithm obtained experimentally. Figure 4 also shows dependences of RMS error of measurement of the UAV position, which corresponds to the lower bound of Cramer-Rao (curve 4). Accumulation of measurements of a trajectory filtration allows to reduce RMS error of location definition of the UAV in comparison with the lower bound of Cramer-Rao (Tovkach and Zhuk 2017b) by 2 to 4 times.

Figure 4
RMS errors of the UAV position estimate using adaptive filter.

Algorithms of estimating the UAV movement parameters using two Kalman filters using models of almost uniform motion j = 2 and motion with maneuver j = 3 were also investigated. Estimation errors of the UAV position of the Kalman filter based on the model j = 2 (Fig. 5) contain systematic components, which in some sections of the trajectory are 2 times higher than the RMS error of the position estimation due to the presence of maneuvers. The developed adaptive algorithm allows to increase the accuracy of the UAV movement parameters estimate on hover and motion without maneuver sections in comparison with the Kalman filter based on the model j = 3 (Fig. 6) more than 2 to 3 times.

Figure 5
RMS errors of the UAV position estimate using Kalman filter based on the model of j = 2.

Figure 6
RMS errors of the UAV position estimate using Kalman filter based on the model of j = 3.

CONCLUSIONS

The optimal algorithm of adaptive filtration of UAV motion parameters based on TDOA measurements (Eqs. 8 to 11) obtained on the basis of a mathematical apparatus of mixed Markov processes in discrete time allows to recurrently calculate the joint PD of a continuously-valued vector and the switching variable at each step. The optimal device that realizes the algorithm (Eq. 8 to 11) has M channels. There is no growth of channels number due to the Markov property of the mixed process.

In the quasi-optimal algorithm of adaptive filtration of UAV movement parameters based on TDOA measurements (Eqs. 8, 11 to 16) obtained by linearization of the equation of measurement (Eq. 3), Gaussian approximation of conditional a posteriori distributions of a continuous-valued component is carried out. At the same time, it estimates the measurement error of the reference sensor, which allows to realize serial processing of measurements and to avoid the inversion of large-dimensional matrices.

Application of a trajectory filtration allows to reduce RMS of errors of definition of the UAV location in comparison with RMS of errors of location definition by TDOA method (the lower bound of Cramer-Rao) by 2 to 4 times. For the reviewed example the developed adaptive algorithm allows to increase the accuracy of estimating the UAV movement parameters on hover and motion without maneuver sections in comparison with the Kalman filter based on the model j = 3 more than 2 to 3 times and to avoid the occurrence of systematic error estimates. At the same time, the probability of recognizing these areas is close to unity.

• FUNDING
  There are no funders to report for this submission.

REFERENCES

• Card BA (2018) How the commercial unmanned aerial vehicle revolution threatens the US threshold. Air & Space Power Journal Spring 2018.
• DroneBouncer (2018) Zashchita ot dronov - nebesnyy rubezh vashey bezopasnosti. DroneBouncer; [accessed 2019 May 24]. http://dronebouncer.com/
  » http://dronebouncer.com/
• Evlanov PA, Zhuk SY (1990) Integration of meters with failures. Radioelectronics and Communications Systems 33(7):49-53. https://doi.org/10.3103/S073527271990070111
  » https://doi.org/10.3103/S073527271990070111
• Friese L, Jenzen-Jones NR, Smallwood M (2016) Emerging unmanned threats: the use of commercially-available UAVs by armed non-state actors. Churchlands Armament Research Services.
• Gemayel NE, Jäkel H, Jondral FK (2014) Error analysis of a low cost TDoA sensor network. Presented at: IEEE/ION Position, Location and Navigation Symposium; Monterey, USA. https://doi.org/10.1109/PLANS.2014.6851484
  » https://doi.org/10.1109/PLANS.2014.6851484
• Makki A, Siddig A, Saad M, Cavallaro JR, Bleakley CJ (2016) Indoor localization using 802.11 time differences of arrival. IEEE Transactions on Instrumentation and Measurement 65(3):614-623. https://doi.org/10.1109/TIM.2015.2506239
  » https://doi.org/10.1109/TIM.2015.2506239
• Nonami K, Kendoul F, Suzuki S, Wang W, Nakazawa D (2010) Autonomous flying robots unmanned aerial vehicles and micro aerial vehicles. Tokio: Springer.
• Rullan-Lara JL, Sanahuja G, Lozano R, Salaza S, Garcia-Hernandez R, Ruz-Hernandez JA (2013) Indoor localization of a quadrotor based on WSN: a real-time application. International Journal of Advanced Robotic Systems 10(1):48:2013. https://doi.org/10.5772/53748
  » https://doi.org/10.5772/53748
• Sathyamoorthy D (2015) A Review of Security Threats of Unmanned Aerial Vehicles and Mitigation Steps. Journal of Defense and Security 6(1):81.97.
• Solodov A, Williams A, Hanaei SA, Goddard B (2018) Analyzing the threat of unmanned aerial vehicles (UAV) to nuclear facilities. Security Journal 31(1):305-324. https://doi.org/10.1057/s41284-017-0102-5
  » https://doi.org/10.1057/s41284-017-0102-5
• Tovkach IO, Zhuk SYa (2017a) Adaptive filtration of parameters of the UAV movement on data from its location calculated on the basis the time difference of arrival method. Presented at: IEEE First Ukraine Conference on Electrical and Computer Engineering; Kiev, Ukraine. https://doi.org/10.1109/ukrcon.2017.8100466
  » https://doi.org/10.1109/ukrcon.2017.8100466
• Tovkach IO, Zhuk SYa (2017b) Recurrent algorithm for TDOA localization in sensor networks. J Aerosp Technol Manag 9(4):489-494. https://doi.org/10.5028/jatm.v9i4.727
  » https://doi.org/10.5028/jatm.v9i4.727
• Tovkach IO, Zhuk SYa (2017c) Adaptive filtration of radio source movement parameters with complex use of sensor network data based on TDOA and RSS methods. Radioelectronics and Communications Systems, 60(12):528-537. https://doi.org/10.3103/s0735272717120020
  » https://doi.org/10.3103/s0735272717120020
• Wallace RJ, Loffi JM (2015) Examining unmanned aerial system threats & defenses: a conceptual analysis. International Journal of Aviation, Aeronautics, and Aerospace 2(4). https://doi.org/10.15394/ijaaa.2015.1084
  » https://doi.org/10.15394/ijaaa.2015.1084
• Wan P, Li Z, Hao B (2016) Time delay estimation of co-frequency signals in TDOA localization based on WSN. Presented at: International Conference on Computer, Information and Telecommunication Systems; Kunming, China. https://doi.org/10.1109/CITS.2016.7546399
  » https://doi.org/10.1109/CITS.2016.7546399
• Zhuk SY (1988) Joint filtering of mixed Markov processes in discrete time. Radioelectronics and Communications Systems 31(1):29-35. https://doi.org/10.3103/S073527271988010066
  » https://doi.org/10.3103/S073527271988010066
• Zhuk SY (1989) Synthesis of digital detector- meters for mixed Markovian processes. Radioelectronics and Communications Systems 32(11):29-34. https://doi.org/10.3103/S073527271989110063
  » https://doi.org/10.3103/S073527271989110063

Edited by

Section Editor: Alison Moraes

Publication Dates

• Publication in this collection
  26 Aug 2019
• Date of issue
  2019

History

• Received
  12 Jan 2018
• Accepted
  21 Nov 2018