Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 4, pp. 915-926, Warsaw 2013

RESULTS OF A RESEARCH PREDICTING THE POSITION OF

AN AIRCRAFT DURING APPROACH AND LANDING USING

THE BESSEL FUNCTION

Marek Grzegorzewski

Air Force Academy in Deblin, Poland

e-mail: marekgrzegorzewski@wp.pl

The article presents theoretical foundations of the position prediction functionality of aGPS
receiver to be used by a pilot in the event of an instantaneous lack of position data occurring
due to various reasons.

Key words: satellite technology, aviation, predication position

1. Navigating an aircraft in space by means of a position potential

Thenumber and character of factors which exert an influence on obtaining the complete naviga-
tion information result in the fact that observation data do not always constitute the complete
information regarding thepreciseassessmentof the locationof agivenaircraft.Evenwhenthe sa-
tellite system is fully available and its configuration proper, some factors (e.g. abrupt changes of
flight parameters ormaneuvers of the aircraft)may prevent its utilization. The above-mentioned
factors necessitate working out additional methods which supplement the process of navigating
an aircraft.Mathematicalmethods are one of themost frequentways of finding solutions to such
problems. In the presentedwork, an alternative filter has been developed which relies mostly on
the following parameters related to motion of an aircraft: position, velocity and position error
statistics. In reality, a cluster of observed position fixes contains all vital kinematic information
of the aircraft.The assumptions of the programuse the concept of the “position potential” based
on the data received from the navigation services. According toNewton’s law of universal gravi-
tation, a freemass particle is attracted by another mass. The acceleration of a particle that has
a givenmass is directly proportional to the gravitation constant G and inversely proportional to
the square of the distance from the attracted mass. Per analogiam, let us regard the statistical
confidence regions (position error ellipsoids) of themost probable position fix as a “source of for-
ce”, which should “attract” a trajectory of a particle passing through them. The potential field
of such a single particle, which in our case is an aircraft, should reflect the observed position fix
and the required force to be exerted on the particle, which shouldmonotonically decrease when
the particle approaches the assumed position fix. When the particle enters the position error
ellipsoid, it will be attracted with a force whosemagnitude will be proportional to the intensity
of the potential. The potential monotonically decreases with the decreasing distance between
the particle and the assumed fix. Furthermore, in order to allow the particle to continuemoving
after the position fix appears, the potential of the attracting source will be dissipated exponen-
tially in time. If we select the position density function (“position potential”), which contains
the dissipation exponent α and the parameter G (positioning uncertainty, corresponding to
Newton’s gravitation constant), we assume that the trajectory of the particle will represent the
real trajectory of the vehicle.

The first attempt at developing a navigation filter that was based on the above assumptions
was made by Inzinga and Vanicek (1985). Their study assumes that the force of the position



916 M. Grzegorzewski

potential field affecting the particle is connected with the probability which describes the fact
that the position fix is included in a two-dimensional error ellipse. Since the equation of motion
of the particle is extremely difficult to be treated analytically, the following model of operation
was selected:
• selecting the function of the position potential,

• formulation of temporal dissipation of the potential field,

• describing the equation of motion,

• solving the equation of motion,

• determining parameters α and G,

• finding the final solution to the navigational problem.

First of all, we have to select the proper potential function. Using the selected function of
position density, we can establish the time-related position potential field for a sequence of
position fixes. Subsequently, we set up amodel of motion of an individual particle and find the
solution to its equation of motion. In order to reflect the changing navigational environment,
the potential function contains certain variable parameters α and G. When the parameters α
and G as well as the initial conditions are not known, the number of possible positions of the
particle is infinite. Because of that, further data is necessary in order to determine the proper
trajectory of the particle. Therefore, we have to set up a “self-learning” procedure for the filter
so that it will determine parameters and conditions based on earlier observations. Parameters
and initialmotion conditions in themotionmodel are related to previous observations bymeans
of “equations by motion”. At this stage, we solve the equation by motion with the assumption
that the initial conditions are known. Then, we have to determine the parameters α and G
and in order to optimize them we will use the least-squares method. The estimation of a given
future position of the particle (aircraft), i.e. making a prediction, is possible owing to themodel
of motion having the determined parameters and on the basis of the present position of the
particle. Error estimation is possible at the verymoment of obtaining the data concerning a new
positionfix, byusing thevalues of the differences between the predicted and the actual positions,
i.e. checking whether the new position fix is located within the position error ellipsoid.

1.1. Navigation model

Weassume that the density function (Plucińska andPluciński, 2000) for a three-dimensional
vector of random position is

Φr0 =
1

K
exp
[

−
1

2
(r−r0)

T
C
−1(r−r0)

]

(1.1)

where r= [x,y,z]T is the particle position vector at the present instant t, r0 = [x0,y0,z0]
T –

vector of particle position at the instant set t0, x,y,z are coordinates of the position vector of
the particle (the aircraft), C – symmetrical covariance matrix (Smirnov, 1967)

C=







D2X cov(X,Y ) cov(X,Z)
cov(Y,X) D2Y cov(Y,Z)
cov(Z,X) cov(Z,Y ) D2Z






(1.2)

and

K =
√

(2π)3detC (1.3)

Let us assume that the particle potential Ui is the basis for the determination of its position
fix

Ui(t)=G(r−r0i)
T
C
−1
i
(r−r0i)e

−α(t−ti) (1.4)



Results of a research predicting the position of an aircraft ... 917

(r−r0i)
TC
−1
i
(r−r0i) is the quadratic form of the error ellipsoid at the instant ti, time t­ ti.

Thepoints located on the error ellipsoid have identical probability distributiondensity. Ci is the
positive covariancematrix of the i-thparticle,whereaspositive parameters αand G (their values
are to be determined) represent, respectively; α – dissipation parameter, with the increasing
value of α parameter, the potential of the particle decreases exponentially; G – uncertainty
of the position fix determination. It is a counterpart of the gravitational constant in Newton’s
attraction.

The potential of the particle is directly proportional to G parameter.

In relation to the potential determined in (1.4), the time-varying t potential created by n
position fixes of a particular particle is expressed as

U =
n
∑

i=1

Ui =Ge
−αt

n
∑

i=1

(r−r0i)
T
C
−1
i
(r−r0i)e

−αti (1.5)

where t­ tn.

It should be noted that the filter (potential) keeps pace with the kinematics of the particle
(aircraft). The time-varying potential field is constantly updated by each newposition fix of the
particle in order to incorporate the most recent information as soon as possible.

Local disturbances (e.g. wind velocity, weather conditions) do not significantly influence the
principal directions of the error ellipsoid. Therefore, we may assume that the principal axes of
the error ellipsoid are parallel to the coordinate axes (X,Y,Z). Consequently, randomvariables
X, Y ,Z are independent and uncorrelated, and the covariance matrix is diagonal

cov(X,Z)= cov(Y,Z)= cov(X,Y )= 0 (1.6)

Equation of motion (1.4), after adopting the following notation

A=2
n
∑

i=1

e−αtiC−1
i

B=2
n
∑

i=1

e−αtiC−1
i
r0i (1.7)

whereA – matrix (3×3), B – vector, is expressed in the coordinates as follows (t­ tn)

ẍ(t)=−GAx
(

x−
Bx

Ax

)

e−αt ÿ(t)=−GAy
(

y−
By

Ay

)

e−αt

z̈(t)=−GAz
(

z−
Bz

Az

)

e−αt
(1.8)

which is a system of inhomogeneous linear differential Bessel equations of the second order (see
Kącki and Siewierski, 2002).

The solution to system of equations (1.8) will be the following

x(t)=
Bx

Ax
+a1J0

(

2

α
e−
α

2
t
√

GAx

)

+a2N0

(

2

α
e−
α

2
t
√

GAx

)

y(t)=
Bx

Ax
+ b1J0

(

2

α
e−
α

2
t
√

GAy

)

+ b2N0

(

2

α
e−
α

2
t
√

GAy

)

z(t) =
Bx

Ax
+ c1J0

(

2

α
e−
α

2
t
√

GAz

)

+ c2N0

(

2

α
e−
α

2
t
√

GAz

)

(1.9)

where the first element of the solution is to be regarded as a particular solution to an inhomo-
geneous differential equation, and both remaining elements, containing the constants, as general
solutions to a homogeneous differential equation.



918 M. Grzegorzewski

In order to determine the integration constants a1, a2, b1, b2, c1, c2, the initial conditions
have to be used, i.e. theymust be added to systemof equations (1.8) (formulation of theCauchy
problem).
Equations (1.9) and (1.10) with known (to be determined) parameters G and α are basic

equations of our study see Grzegorzewski (2011).
Using the initial condition, for tn =0

r(0)=







xn
yn
zn






ṙ(0)=







ẋn
ẏn
żn






(1.10)

we compute the constants a1, a2, b1, b2, c1, c2

a1 =−
π

α

[(

xn−
Bx

Ax

)

√

GAxN1

(2

α

√

GAx

)

− ẋnN0

(2

α

√

GAx

)]

a2 =
π

α

[(

xn−
Bx

Ax

)

√

GAxJ1

(2

α

√

GAx

)

− ẋnJ0

(2

α

√

GAx

)]

b1 =−
π

α

[(

yn−
By

Ay

)
√

GAyN1

(2

α

√

GAy

)

− ẋnN0

(2

α

√

GAy

)]

b2 =
π

α

[(

yn−
By

Ay

)
√

GAyJ1

(2

α

√

GAy

)

− ẋnJ0

(2

α

√

GAy

)]

c1 =−
π

α

[(

zn−
Bz

Az

)

√

GAzN1

(2

α

√

GAz

)

− ẋnN0

(2

α

√

GAz

)]

c2 =
π

α

[(

zn−
Bz

Az

)

√

GAzJ1

(2

α

√

GAz

)

− ẋnJ0

(2

α

√

GAz

)]

(1.11)

Because of initial condition (1.10), the computed integration constants a1, a2, b1, b2, c1, c2
represented by (1.11) exist within the time interval [tn, tn+1] and [x,y,z] – vector of the particle
(aircraft) position, [m]; [ẋ, ẏ, ż] – vector of the particle (aircraft) velocity, [ms−1].

1.2. Optimization of α and G parameters

In order to optimize α and G parameters, let us analyze the following problem, consisting
in solving the equation of motion (T ­ tn) (Grzegorzewski, 2011)

ẍ(t)=−G(Axx−Bx)e
−αt ÿ(t)=−G(Ayy−By)e

−αt

z̈(t)=−G(Azz−Bz)e
−αt

(1.12)

where the vector Ar is represented by

Ax =2
n
∑

i=1

eα(ti−tn)pxi Ay =2
n
∑

i=1

eα(ti−tn)pyi Az =2
n
∑

i=1

eα(ti−tn)pzi (1.13)

and the coordinates of the vector B are represented as

Bx =2
n
∑

i=1

eα(ti−tn)pxix0i By =2
n
∑

i=1

eα(ti−tn)pyiy0i

Bz =2
n
∑

i=1

eα(ti−tn)pziz0i

(1.14)

with the initial conditions set for tn =0 as in Eq. (1.10).
The vector of the particle position is described by Eqs. (1.9) (Grzegorzewski, 2011).



Results of a research predicting the position of an aircraft ... 919

Our objective is to optimize the parameters α and G so that the function of two variables

f(α,G) =
n
∑

i=1

[r(ti)−r0i]
2 =

n
∑

i=1

{

[x(ti)−x0i]
2+[y(ti)−y0i]

2+[z(ti)−z0i]
2
}

=min (1.15)

where: x0i,y0i,z0i is the position of the particle (aircraft) at the instance ti (data obtained from
a satellite), x(ti),y(ti),z(ti) – position of the particle (aircraft) at the instance ti computed
from (1.9).
By substituting the optimized parameters α and G into equations (1.9) and (Grzegorzewski,

2011) we will obtain from themathematical model the vector of the particle (aircraft) position
and the vector of its velocity at the instant t.

2. Comparing various methods of predicting the aircraft position

Our task was to create an algorithm, which may be implemented and then used in practice, in
order to predict subsequent positions of an aircraft.
The results of our work are the following:

• Three different algorithms, written as MS Excel macros, which predict the position of an
aircraft at the subsequent instances on the basis of its positions at the previous instances.

• Comparison of the efficiency of the created algorithms when using various values of the
parameters.

• Suggestions concerning the proper selection of an algorithm and its implementation.

2.1. Methods under comparison

The starting point was the method described by Grzegorzewski (2011). Its main advantage
is the fact that it has been created on the basis of real physical assumptions. This method is
henceforth called αGnmethod – from its two key parameters. In our studies, implementation of
other methods was also taken into consideration. Having analyzed their theoretical correctness,
we carried out their preliminary tests on the basis of the real data contained in the MS Excel
sheet, data.xls. Two methods proved to be worth further testing. They represent two different
approaches:

a) Establishing the relationship between the coordinates of the aircraft position fix and time.

b) Establishing the relationship between the aircraft position at a given instant of time and
its position at the previous time instances.

In the next part of this reportwe provide a brief discussion of the threemethods and the results
of the tests conducted in order to assess their usefulness.

2.1.1. αGn method

Technically, its essence lies in expressing the solution of an ordinary second order linear
differential equation as a function of two parameters α and G, and, subsequently, selecting the
optimum values of those parameters which will minimize the error function F0. Knowing those
values facilitates the determination of thepredicted values.The solution of the equation contains
Bessel special functions.
Function f(α,G) represents the sum of squares of the distances between the theoretical po-

ints and the observed points at n previous time instances. Initial terms of the series expansion
of the Bessel functionwere used for the approximation of that function. Due to high complexity
of the formulae describing the function f(α,G), applying traditional optimization methods was



920 M. Grzegorzewski

not possible. Therefore, we decided to use for its optimization a certain version of an evolu-
tionary algorithm developed by ourselves. The implementation of that algorithm is, however,
time-consuming, and increasing the number of previous time instances n causes a considerable
increase of duration of the algorithm operation. The parameters affecting the course of action
of the evolutionary algorithm are: population size, number of discarded candidate solutions
(individuals), width (step) of the local method, number of steps, and finally ranges: α,G.

2.1.2. lin n method

According to our experience and the data that we obtained, predicting the height h is of
crucial importance whereas predicting the remaining coordinates X and Y is considerably less
important. Therefore, both this method and the next are being discussed in the context of
predicting the height, but the discussionmay, if necessary, be repeated separately for X and Y .
This is one of the simplest methods. It assumes a linear relationship between the height and
time

h(t) = b0+ bit (2.1)

The coefficients b0, bi are determined with the least squares method on the basis of the course
of the flight so far (n measurements, n> 2), and they are used for making a prediction. After
preliminary tests, other methods of determining a functional relationship between the height
and time were discarded as being substantially less effective.

2.1.3. baryc n method

Thismethod is in fact one of the versions of the ARIMA(p,d,q) method, namely A(1,1,0),
and its description is expressed as: hi =β1hi−1+β2hi−2 with β1+β2 =1.

The coefficients β1, β2 are determined with the least squares method on the basis of the
course of the flight so far (n measurements, n > 2) and used for making a prediction. After
preliminary tests, other methods representing the ARIMA type were discarded as being sub-
stantially less effective. For similar reasons, two other methods of a similar type were discarded
after preliminary tests.

2.1.4. In-flight test results

Results representative for the data set included inMSExcel sheet, data.xls (flight of theCes-
sna aircraft) are presented in Fig. 1a. The data was used for comparing themethodsmentioned
above and for optimizing the α and G parameters.

Fig. 1. (a) Height in the data.xls sheet; (b) height in the 4 minutes 92.xls sheet



Results of a research predicting the position of an aircraft ... 921

Table 1 contains the results for the data set fromMSExcel sheet, 4minutes 92.xls (Fig. 1b).
Predictions were examined from 10th to 240th second, i.e. at 231 instants. Similar calculations
were performed for data.xls (Fig. 1a).

Table 1.Parameters alpha and G determined for αG2 and αG3methods

α G

αG2method

Mean 13.01 0.97

Stand. deviation 5.78 1.00

Minimum 1.00 0.010003

Maximum 19.97 2.993

αG2method

Mean 14.22 0.12

Stand. deviation 5.25 0.22

Minimum 1.19 0.010003

Maximum 19.98 1.95

Numerical experiments were performed in the actual flight, shown in Fig. 1, using the follo-
wingmethods:αGn, lin n, and baryc n.The last twomeasurementswere taken into consideration
(in the case of αGn method, this means that the memory of the potential was reduced to 2 se-
conds), and the predicted position of the aircraft was calculated for 2 or 3 seconds forward (due
to the complexity of calculations –Bessel functions). At the same time, particularmethodswere
compared and analyzed for the number of times that a givenmethodwas the best. On the basis
of the numerical experiments, it is believed that the followingmethods should be used in further
in-flight tests: αG2, αG3, and lin 2. In the event of substantial flight disturbances, αG2 and
αG3 methods “adapts to the actual flight trajectory” faster.

3. Position prediction – results

Experiments concerning the position prediction were conducted during flights designated as
Flight 12, and Flight 22.

3.1. Flight 12

Flight 12was performed on 12 September 2011, between 9:02:40s and 9:04:00s.On the basis
of the SV arrangement, PDOP=2.25. Waypoints of the glide path which is shown in Fig. 1a
were determined once a second. The following points were planned for the analysis:

• reference point 3,

• reference point 6,

• reference point 9.

Figure 2 shows the parameters of the glide path inFlight 12. The aimof the graph is to show
the trajectory of theglide pathof anaircraft (coordinates recordedonce every 10 seconds). In the
three points thatwere planned to be shown in the graph, theGPS system stopsworking and the
positionprediction system takes over. The comparison of the positionprediction result at point 3
(3rd second of the aircraft flight on the glide path) with the position obtained from the GPS
shows that the real position of the aircraft was 5 meters higher than the position computed by
means of prediction. The directional deviation from the glide path was also recorded – 7meters
to the right.



922 M. Grzegorzewski

Fig. 2. Glide path graph – Flight 12

After the subsequent second of flight, the GPS is deactivated again. At point 6 (6th second
of flight), the comparison of the position prediction result with the position obtained from the
GPS shows a 7-metre height difference between the real position of the aircraft and the position
computed by means of prediction. The position prediction system maintains the tendency to
compute the parameters of the subsequent positions above the real height. Taking into conside-
ration the safety of flight and landing, such a tendency guarantees a safe landing.

Table 2.Geographical coordinates – Flight 12

GPS
coordinates

Air- GPS Point coordinates Point coordi-
Point speed hea- – prediction with A nates – linear B

[km/h] ding the heading prediction (′′)

4 51◦34′25.607′′N, 227 263◦

21◦53′51.036′′E

5 51◦34′22.008′′N, 228 256◦ 51◦34′23.159′′N, 54m 51◦34′26.543′′N, 156m
21◦53′21.443′′E 21◦53′20.508′′E 21◦53′17.987′′E

6 51◦34′14.735′′N, 226 245◦ 51◦34′16.68′′N, 72m 51◦34′23.087′′N, 320m
21◦52′51.528′′E 21◦52′22.979′′E 21◦52′42.672′′E

7 51◦34′0.984′′N, 224 190◦

21◦52′49.979′′E

A – Prediction with the heading – value of deviation (′)
B – Linear prediction value of deviation (′′)

Table 2 contains the coordinates of waypoints in the aircraft trajectory in the turn into
the landing heading. Control points 5 and 6 were selected for research purposes. The point
coordinates are computed every 10 seconds of flight. Analysis of the flight during the turn into
the landing heading is particularly significant from the point of view of ensuring the safety of
the approach to landing at the destination aerodrome. Correct execution of this flight phase in
the IFR conditions guarantees safe landing of the aircraft. The commencement of that flight
phase at the exactly prescribed point in the air and maintaining flight parameters during the
turn, i.e. airspeed, bank angle and the turn rate, have a considerable influence on the final result
of executing the turn.

Figure 3 illustrates a segment of the trajectory during the turn into the landing heading.
At point 4, the GPS system is deactivated and the position prediction system takes over. After



Results of a research predicting the position of an aircraft ... 923

10 seconds of flight to point 5, the difference between the real position of the aircraft and its
position obtained from “linear prediction” (black color, point 5′′) is 155 meters. In the system
that hasmemorized the tendencies of heading changes during the flight and has taken into acco-
unt those changes in the position predictionwith the heading, the difference between position 5′

(blue color) and the real position (point 5) is only 40meters. Continuation of the flightmainta-
ining the same parameters results in further deviation between the aircraft computed position
from its real trajectory. The distance between point 6′′ and point 6 is 309 meters whereas the
distance between points 6′ and 6 is 67 meters. On the basis of the assessment of the position
prediction system efficiency in the turn into the landing heading in IMC (Instrument Mete-
orological Conditions) flights, a conclusion should be made that “prediction with the heading”
facilitates establishing the aircraft into the runway heading.

Fig. 3. Flight trajectory during the turn into landing heading – Flight 12

3.2. Flight 22

Flight 22 was performed on 22 September 2011, between13:15:05s and 13:16:45s. On the
basis of the SV arrangement, PDOP=2.1. The following points were planned:

• reference point 2,

• reference point 4,

• reference points 7 and 8,

• reference point 10.

The parameters of the glide path in Flight 22 are shown in Fig. 4. The aim of the graph is
to show the trajectory of the glide path of an aircraft (coordinates recorded once a second) and
present the results of the position prediction system. In the four points that were planned to
be shown in the graph, theGPS system stops working and the position prediction system takes
over. The comparison of the position prediction result at point 4 (4th second of the aircraft
flight on the glide path) with the position obtained from the GPS shows that the real position
of the aircraft was identical with the position computed bymeans of prediction. the directional
deviation from the glide path was also recorded – 2 meters to the right of the runway heading.
In the 6th second of flight, the GPS is deactivated again. At point 8 (8th second of flight),
the height difference between the position obtained by use of the position prediction and the
position obtained from the GPS was−2 meters.
Figure 5 illustrates a segment of the trajectory during the turn into the landing heading.

At point 1, the GPS system is deactivated and the position prediction system takes over. After
10 seconds of flight to point 2, the difference between the real position of the aircraft and its



924 M. Grzegorzewski

Fig. 4. Glide path graph – Flight 22

Table 3.Geographical coordinates – Flight 22

GPS
coordinates

Air- GPS Point coordinates Point coordi-
Point speed hea- – prediction with A nates – linear B

[km/h] ding the heading prediction (′′)

1 51◦32′26.735′′N, 222 267◦

21◦51′36.144′′E

2 51◦32′25.583′′N, 223 275◦ 51◦32′25.367′′N, 7m 51◦32′24.719′′N, 27m
21◦50′58.271′′E 21◦50′58.163′′E 21◦50′58.127′′E

3 51◦32′28.608′′N, 221 302◦

21◦50′20.651′′E

4 51◦32′45.599′′N, 224 28◦ 51◦32′43.007′′N, 178m 51◦32′36.383′′N, 393m
21◦49′53.436′′E 21◦49′45.084′′E 21◦49′39.215′′E

5 51◦33′8.208′′N, 225 55◦

21◦50′4.956′′E

6 51◦33′24.335′′N, 225 74◦

21◦50′27.743′′E

7 51◦33′34.848′′N, 223 91◦ 51◦33′36.72′′N, 59m 51◦33′40.319′′N, 188m
21◦51′4.571′′E 21◦51′4.14′′E 21◦51′0.36′′E

8 51◦33′34.56′′N, 220 100◦ 51◦33′40.319′′N, 179m 51◦33′57.167′′N, 710m
21◦51′45.143′′E 21◦51′44.495′′E 21◦51′39.276′′E

9 51◦33′57.167′′N, 190 107◦

21◦51′39.276′′E

10 51◦33′17.856′′N, 180 106◦ 51◦33′19.115′′N, 40m 51◦33′20.159′′N, 74m
21◦52′52.284′′E 21◦52′52.608′′E 21◦52′53.327′′E

A – Prediction with the heading – value of deviation (′)
B – Linear prediction value of deviation (′′)

position obtained from“linear prediction” (black color, point 2”) is 27meters. In the system that
has memorized the tendencies of heading changes during the flight and has taken into account
those changes in thepositionpredictionwith theheading, thedifferencebetweenposition 2′ (blue
color) and the real position (point 2) is only 7 meters. After that the flight is GPS-controlled
again. At point 3, theGPS system is deactivated and the position prediction system takes over.
The accuracy of prediction is assessed at point 4. The distance between point 4′′ and point 4 is



Results of a research predicting the position of an aircraft ... 925

393meters whereas the distance between points 4′ and 4 is 178meters. After exiting the turn to
establish on landing heading, point 8′′ is 710meters form point 8 whereas point 8′ is 179meters
from 8.After exiting the turn to establish on landing heading at point 10, point 10′′ is 74meters
form point 10 whereas point 10′ is 40 meters from 10.

Fig. 5. Flight trajectory during the turn into landing heading – Flight 22

4. Summary and conclusions

On the basis of the conducted numerical experiments, we hold the view that themethods to be
implemented in further in-flight tests are:

a) αG2 and αG3,

b) linear 2.

When there are no sudden changes in the course of the flight, the accuracy of the methods
mentioned above is considerably higher than the accuracy of the remaining methods. However,
only after implementing them in a real instrument, a GPS receiver, the real operational time
efficiency of each of themmay be assessed. According to our predictions, the linear 2method is
a little faster than the αG2, which in turn is faster than αG3. In the event of substantial flight
disturbances, none of those methods, or any other method known to us, may be used without
running a great risk of making a serious error. During the flight, making the position prediction
ought to be repeated once every second on the basis of the known position fix of the aircraft
in the previous 2 or 3 seconds. The prediction may be used as necessary (lack of GPS data)
provided that there are no sudden changes in the flight parameters. In the event of such changes,
the prediction should not be utilized.
The accuracy of the prediction decreases with time. Therefore, the pilot ought to be warned

after a given time (we recommend 3 seconds) that the position of the aircraft being shown
may differ considerably from its real position. Glide paths in Figs. 2 and 4 show the differences
between real and predicted positions of the aircraft in the vertical plane. The differences in the
horizontal plane are illustrated in Figs. 3 and 5, allowing the determination of themagnitude of
the corrections necessary tobeput in theflightpath inorder to obtain theheadingcorresponding
to the runway axis.
Glide paths inFigs. 2 and 4 showthat themagnitude of corrections tobemade to the aircraft

altitude, meets the requirements of Category I and II (according to ICAO – Requirements for
Navigation Systems in the Field of Reliability, Availability and Accuracy).



926 M. Grzegorzewski

The presented solution facilitates the performance of the aircraft landingwithin Category II
(ICAO). There is still a need, however, to continue research on the use of position prediction.

References

1. Grzegorzewski M., 2005, Navigating an aircraft by means of a position potential in three-
dimensional space,Annual of Navigation, 9, ISSN 1641-9723

2. GrzegorzewskiM., 2011,Matematycznymodel nawigacji statku powietrznego z użyciem wektora
położenia i prędkości oraz jego zastosowanie w odbiorniku GPS do predykcji pozycji w trójwymia-

rowym układzie współrzędnych, Dęblin, ISBN 978-83-60908-65-5

3. Inzinga T., Vaniček P., 1985, A two-dimensional navigation algorithm using a probabilistic
force field, Proceedings of Third International Symposium on Inertial Technology for Surveying
and Geodesy, Banff, Canada

4. Kącki E., Siewierski L., 2002,Wybrane działymatematyki wyższej z ćwiczeniami, Łódź,Wyższa
Szkoła Informatyki

5. Kościelniak P., Ombach J., 2010, Sprawozdanie – zadanie badawcze nr 1 (Projekt
0028/B/T00/2008/35), Biblioteka Główna Wyższa Szkoła Oficerska Sił Powietrznych U-3941/1,
Dęblin

6. MatwiejewN.M., 1972,Metody całkowania równań różniczkowych zwyczajnych,Warszawa,PWN

7. OmbachJ., 2004,Somealgorithmsofglobaloptimizations,Post-ConferenceMaterials Jagiellonian
University, Cracow

8. Platt C., 1981, Problemy rachunku prawdopodobieństwa i statystyki matematycznej, Warsaw,
PWN

9. Plucińska A., Pluciński E., 2000,Probabilistyka. Macierz kowariancji, Warszawa,WNT

10. Smirnow W.I., 1967,Matematyka Wyższa, Warszawa, PWN

Manuscript received January 6, 2013; accepted for print March 12, 2013