IntechOpen

Home > Books > Introduction and Implementations of the Kalman Filter

Open access peer-reviewed chapter

Tuning of the Kalman Filter Using Constant Gains

Written By

Mudambi R. Ananthasayanam

Submitted: 05 March 2018 Reviewed: 02 October 2018 Published: 23 November 2018

DOI: 10.5772/intechopen.81795

Introduction and Implementations of the Kalman Filter IntechOpen
Introduction and Implementations of the Kalman Filter Edited by Felix Govaers

From the Edited Volume

Introduction and Implementations of the Kalman Filter

Edited by Felix Govaers

Book Details Order Print

Chapter metrics overview

2,077 Chapter Downloads

View Full Metrics
Advertisement

Abstract

For designing an optimal Kalman filter, it is necessary to specify the statistics, namely the initial state, its covariance and the process and measurement noise covariances. These can be chosen by minimising some suitable cost function J . This has been very difficult till recently when a near optimal Recurrence Reference Recipe (RRR) was proposed without any optimisation but only filtering. In many filter applications after the initial transients, the gain matrix K tends to a constant during the steady state, which points to design the filter based on constant gains alone. Such a constant gain Kalman filter (CGKF) can be designed by minimising any suitable cost function. Since there are no covariances in CGKF, only the state equations need to be propagated and updated at a measurement, thus enormously reducing the computational load. Though CGKF results may not be too close to those of RRR, they are acceptable. It accepts extremely simple models and the gains are robust in handling similar scenarios. In this chapter, we provide examples of applying the CGKF by ancient Indian astronomers, parameter estimation of spring, mass and damper system, airplane real flight test data, ballistic rocket, re-entry of space object and the evolution of space debris.

Keywords

  • adaptive EKF
  • reference recursive recipe
  • maximum likelihood
  • Cramer Rao bound
  • constant gain Kalman filter

1. Introduction to Kalman filter

The simplest formulation of a Kalman filter [1] is when the state and measurement equations are both linear. However, Kalman filter has found its greatest application for non-linear systems. A typical continuous state with discrete measurements in time forming a non-linear filtering problem can be written as

x k = f x k − 1 Θ u k − 1 + w k , E1
Z k = h ( x k Θ ) + v k , k = 1 , 2 , 3 , … , N E2

where ‘x’ and ‘Z’ are, respectively, the state and measurement equations of size (n × 1) and (m × 1); u is the control input and Θ the parameter vector of size (p × 1) and ‘f’ and ‘h’ are non-linear functions. The process ‘w’ and measurement ‘v’ noises are respectively of size (n × 1) and (m × 1). These are assumed to be zero mean with covariances Q and R and their sequences are uncorrelated with each other. The states may not be in general observable but the measurements should be related to the states. In many applications for linear systems, if the unknown parameters Θ are treated as additional states, then the linear system of equations becomes non-linear. In such cases, the extended Kalman filter (EKF) formulation can be written as

x k Θ k = f x k − 1 Θ k − 1 u k − 1 Θ k − 1 + w k 0 E3

or equivalently

X k = f X k − 1 + w k E4
Z k = h X k + v k , k = 1 , 2 , … , N E5

where ‘ X ’ is the augmented state of size ((n + p) × 1). The control symbol ‘u’ is not shown for brevity. The formal solution for the above filtering problem can be summarised following Brown and Hwang [2] as

Initial state estimate X 0 0 = X 0 = E X t 0 , E6
Initial state covariance matrix P 0 0 = P 0 = E X 0 – X t 0 X 0 – X t 0 T E7
Prediction step : X k k − 1 = f X k − 1 k − 1 , E8
P k k − 1 = F k − 1 P k − 1 k − 1 F k − 1 T + Q k E9

We assume that X k k − 1 and P k k − 1 denote the estimates of the state and its covariance matrix, respectively, at time index k, based on all information available up to and including time index k−1. Then, we seek to update the state value from X k k − 1 to X k k using the measurement Z k with uncertainty denoted by R k based on the value of K k called the Kalman gain such that the updated covariance P k k having the individual terms along its major diagonal is a minimum, leading to

Update step : K k = P k k − 1 H k T H k P k k − 1 H k T + R k − 1 E10
X k k = X k k − 1 + K k Z k – h X k k − 1 = X k k − 1 + K k ν k E11
P k k = I − K k H k P k k − 1 E12

with P denoting uncertainty, F k − 1 is the state Jacobian matrix (∂f/∂X) evaluated at X = X k − 1 k − 1 . X k k − 1 denotes the estimate at t(k) based on the process dynamics between t(k−1) and t(k) but before using the measurement information. The measurement Jacobian H k = (∂h/∂X) is evaluated at X = X k . The difference between the actual measurement and the predicted model output

ν k = Z k – h X k k − 1 E13

is called the innovation. The importance of the innovation following white Gaussian for filter performance was brought out by Kailath [3]. When the innovation is white, it means all the information has been extracted from the data and no further information is left out, thus both the models and the algorithm have done their best job.

There are thus five basic filter operations, namely: (i) the state propagation, (ii) the covariance propagation, (iii) Kalman gain evaluation, (iv) the state update and (v) the covariance update. The first and fourth refer to sample values and the second, third and fifth refer to the population characteristics. At any given time point, the statistical combination of the two estimates, one from state and the other from measurement equation, are formal if only the covariances denoting their uncertainties are available. Thus the states and the covariances at all times can be estimated if the initial X0 and P0 as well as Q(k) and R(k) are specified over time, but this is not easy. These have to be specified over a time span in order to match and minimise a cost function based on the innovation, or any other in some best possible sense. A well-known criterion is the method of maximum likelihood estimation (MMLE). When Q ≡ 0, the Kalman gain matrix is zero and the technique is called as the output error method (MMLE-OEM). When Q > 0, the method is called as filter error method [4, 5]. For optimal design of the Kalman filter, the innovation follows a white Gaussian distribution which is operationally equivalent to minimising the cost function

J = 1 / N Σ ν k [ H k P k k − 1 H k T + R k ] − 1 ν k T E14
= 1 / N Σ ν k R − 1 ν k T or E15
= J X 0 P 0 Q R Θ or E16
= J X 0 Θ K traded for P 0 Q and R E17

based on summation over all the N measurements. Thus, the filter has to be tuned or in other words should solve for either the statistics X0, Θ, P0, Q and R, in Eq. (16), or for X0, Θ and K in Eq. (17). Of course, there have been many number of cost functions used in the literature, the only constraint being all should lead to reasonable answers that are acceptable. The Q ≡ 0 case leads to an optimisation of the cost function. If Q > 0, then the filter approach becomes compulsive and generally the cost function is forgotten and mostly the filter statistics are tuned manually to obtain the results.

One can see straightaway that the structure of the above cost function J becomes different due to the change of variables (different combination of the statistics can lead to the same gain!); and hence whatever the results are generated, they will be different but have to be within reasonable limits. In the RRR studies [6, 7, 8, 9], many typical cost functions have been stated to bring home the above point. One can be around the true answer but not at the answer which is not known due to the occurrence of the random unknown sequence of the noise distribution in the data. Hence, estimation theory being an inverse problem, the results are subjective and not objective as many claim. In fact, the whole of statistics is subjective from the beginning to the end and thus the results generated can be stretched to any limit but have to be meaningful, acceptable and useful for further use. As is well known, inverse problems do not have unique answers, more so with randomness being introduced. Unless the above sequence of noise distribution can be worked out correctly, there is no way to get the true answers. Thus, the statistical percolation effect affects all the unknowns in any estimation theory. This should be kept in mind to understand any result based on filter statistics or filter gains in Kalman filtering.

1.1 The competence and beauty of the Kalman filter

The earliest Kalman filter formulation by Kalman [1] dealt with state estimation. However, it has grown at present to handle myriad other scenarios such as state and parameter estimation, data fusion and many more. The Kalman filter can ably estimate or account for time-invariant or time-varying (i) unknown, (ii) inaccurately known or (iii) even unmodellable structure of the state and measurement model equations and the parameters in them as also (iv) the deterministic or random inputs and by accounting for them suitably as process and measurement noises. It can compensate even for computational errors during the entire filter operation.

Further, the state and covariance updates at a measurement depend only on the covariances of the state and the measurements and not their probability distribution (!). Hence, after assimilating the measurement information, the update is subtly reset to follow a Gaussian distribution. Thus, the use of only the estimate and covariance all over the filter tacitly implies (one can however improve the numerical values of the estimate and covariance in non-linear problems) the state and measurement variables are all distributed or approximated as quasi-Gaussian. Hence, with all such subjective features, the final result can only be an answer rather than a true or unique answer. All the above have to be checked for the consistency of the whole process of modelling, convergence of the numerical algorithm and other consistency checks among the variables occurring in the filter as discussed in [6, 10, 11, 12].

1.2 Use of filter statistics in designing the Kalman filter

Assuming that the measurements are available at N discrete time instants, the normalised innovation cost function J fundamental to the Kalman filter as suggested by Sorenson [13] is defined as

J = 1 / N Σ ν k R − 1 ν k T E18

where R is the covariance matrix of the innovation H k P k k − 1 H k T + R k . Here, R which is a function of P0, Q and R varies with time. The estimation of the system parameters X 0 , Θ , P 0 , Q   and R is called filter design or filter tuning as mentioned earlier. Though there are many techniques for adaptively tuning the filter statistics [14], the recent RRR [6, 7, 8, 9] or the heuristic approach of Myers and Tapley [15] for Q , and R , and of Gemson [16] and Gemson and Ananthasayanam [17] for P 0 are perhaps the simplest ones.

1.3 Use of filter gains and design of constant gain Kalman filter (CGKF)

In the constant Kalman gain formulation (in discrete form), the update step

K k = P k k − 1 H k T H k P k k − 1 H k T + R k − 1 E19
X k k = X k k − 1 + K k Z k – h X k k − 1 = X k k − 1 + K k ν k E20
P k k = I − K k H k P k k − 1 E21

gets simplified to determine only the constant gain matrix K(k) by subsuming P 0 , Q and R . Hence, there are no covariance equations for propagation at all, thus enormously reducing the numerical effort and time. The constant filter gain approach is less explored than the filter statistics approach. Many attempts have been made by Wilson [18], Cook and Dawson [19], Grimble et al. [20], Kobayashi [21] and Liu et al. [22]; but these are not simple, except the modified gain extended Kalman filter (MGEF) proposed by Song and Speyer [23]. Gelb [24] and Sugiyama [25] consider the CGKF approach but their method of tuning the desired filter gain parameters is manual. The present rational procedure is based on a suitable normalised innovation cost function as in space debris [26, 27, 28, 29] and many other illustrative examples to follow.

When a regular Kalman filter using the filter statistics operates on the data in general, it turns out that after the initial transients the Kalman gain matrix tends to a constant value. Such a feature has been noticed in the tracking of ballistic rockets by Sarkar [30], evolution and expansion of the space debris scenario and prediction of re-entry objects [27–29]; for parameter estimation of dynamic systems by Viswanath [31]; in rendezvous and docking studies [32, 33] and total electron content in the ionosphere [34, 35] and in integration of GPS and INS [36, 37], target tracking in wireless sensor networks by Yadav et al. [38] and many more. Due to limited space, only the first three will be discussed in this chapter. This observation provides a possible approach in which instead of tuning the usual Kalman filter statistics for X 0 , P 0 , Q and R , in general, a smaller number of Kalman gain elements can be worked out. This constant gain matrix K can be obtained by making the above normalised innovation cost function equal to the number of measurements by assuming the above R in Eq. (18) to be a constant. Then, the R can be estimated as if it is the estimation of pure noise only as is the case of MMLE-OEM [4, 5]. There could be some differences in the gain values obtained from the adaptively or manually tuned P 0 , Q and R and the constant gain approach due to the relative periods of transients and the steady-state conditions. Though the results may not be as close to the optimum, the estimates are generally acceptable. The present examples mostly utilise the genetic algorithm [39, 40] to minimise the cost function J and obtain optimum K. However, before applying the constant Kalman gain approach, it is desirable to carry out extensive studies using any adequate adaptive filtering technique such as RRR. A comparison of the results from the adaptive technique and the constant Kalman gain approach provides confidence in the latter approach. The following sections provide example applications of designing constant gain Kalman filters.

Advertisement

2. Ancient Indian astronomers implicitly using the constant Kalman gain approach

Ancient Indian astronomers needed to calculate the position of celestial objects like Sun, Moon and other planets for timing the Vedic rituals. But their predicted positions changed over many centuries due to unmodelled or unmodellable causes. The philosophy of ‘change, capture and correct’ is the one that is followed in the Kalman filter. The ancient Indian astronomers had understood the above philosophy. They used the above concept to update the parameters for predicting the position of celestial objects based on measurements carried out at various time intervals which can be stated as.

Updated parameter = Earlier parameter + Some quantity × Measured − Predicted Position of the celestial object E22

They could not have done it in any other way to update the earlier parameters called ‘cannons’. The ‘some quantity’ as we will see later on is the Kalman gain. There were no frills and fashion of distributions and the spread to infer uncertainties after combining statistically one estimate in certain units with another estimate in another unit but related to the former. The measured longitude of the celestial object is different from the state that is updated, which is the number of revolutions in a yuga just as state and measurements are in general different in many Kalman filter applications!

Billard [41, 42] had stated that if the elements of Aryabhata are now wrong, they must have been accurate when he was living. Then, newer astronomical elements can be established based on the earlier astronomical elements and the new observations of the present time. Billard [41, 42] provides many cannons starting from around AD 500 by Aryabhata to AD 1600 based on later measurements carried out (over many years or even decades!) to make the predicted position of the objects consistent with new observations. One such canon around AD 898 shows a very high accuracy valid over a larger number of centuries. Sarma [43] quotes such revisions over a period of time. Nilakantha (around AD 1443) had stated that the eclipses cited in Siddhanthas as well as those currently observable can be studied and future eclipses can be predicted (extrapolation!). Also, for the eclipses occurring at other longitudes and latitudes, the predictions can be perfected (data fusion!). Based on these, the past eclipses of one’s own place can be refined equivalent to ‘smoothing’! It is strongly urged that research is undertaken on Billard’s work available in French.

Advertisement

3. Typical parameter estimation studies

In order to illustrate the ability of CGKF, we consider the parameter estimation of a spring, mass and damper (SMD) system with a weak non-linear spring constant and also a real flight test data of an airplane.

3.1 Analysis of spring, mass and damper (SMD) system

The SMD system with weak non-linear spring constant in continuous time (t) is governed by the equations

x 1 ̇ t = x 2 t E23
x 2 ̇ t = − Θ 1 x 1 t − Θ 2 x 2 t − Θ 3 x 1 3 t E24

where x 1 and x 2 are the displacement and velocity states. The ‘dot’ represents differentiation with respect to time (t). The unknown parameter vector Θ = [ Θ 1 , Θ 2 , Θ 3 ]T has the true value Θtrue = [4, 0.4, 0.6]. Θ 3 being a weak parameter, it does not affect the system dynamics much and hence its estimation also has more uncertainty. The complete state vector X = [ x 1 , x 2 , Θ 1 , Θ 2 , Θ 3 ]T of size [(n + p) × 1] which in this case is (5 × 1). The measurement equation is given by.

Z k = H X k E25

where H = 1 0 0 0 0 0 1 0 0 0 is the measurement matrix of size (m × (n + p)) where m = n = 2 and p = 3.

At a measurement with the additional terms to assimilate the measurement data, the above equations become

x 1 ̇ t = x 2 t + K 11 Z 1 – x 1 + K 12 Z 2 – x 2 E26
x 2 ̇ t = − Θ 1 x 1 t − Θ 2 x 2 t − Θ 3 x 1 3 t + K 21 Z 1 – x 1 + K 22 Z 2 – x 2 E27

where K 11 , K 12 , K 21 and K 22 are the elements of the constant Kalman gain matrix K. These along with the parameter vector Θ have to be estimated by minimising the earlier mentioned innovation cost function J . A total of N = 100 simulated measurement data are generated with initial state conditions 1 and 0, respectively, in steps of dt = 0.1 s between time 0 and 10 s.

3.2 Remarks on the SMD parameter estimation

The CGKF were based on 25 Newton-Raphson iterations and RRR results were generated based on 100 iterations for each data set (for obtaining generally four-digit accuracy though not necessary) and are compared in Table 1 below. The parameter values are to be read as for [ Θ 1 , Θ 2 , Θ 3 ]. For Q ≡ 0 and Q > 0 cases, the mean and standard deviation of the parameter estimates from CGKF and RRR are by and large close. In fact, the CGKF estimates are generally within about 1σ of the CRB values given by RRR. In the RRR with constant P0, Q and R, the filter is able to follow the system fairly well due to the time-varying gains providing near optimal solution. But in the CGKF approach, since the gains are constant, the filter is unable to follow the system model as well as by RRR. Thus, CGKF follows a slightly different dynamical model than RRR and hence their results are somewhat different.

SMD SYSTEM: For CGKF, the standard deviation (STDV) is based on parameter estimates and for RRR, the Cramer-Rao Bound (CRB ~ STDV) is based on filter covariance averaged over 100 simulations
Results for the three parameters
Case Mean STDV
CGKF Q ≡ 0 3.9982   0.3994   0.5948 0.0377   0.0112   0.0954
RRR Q ≡ 0 4.0028   0.4000   0.5921 0.0242   0.0055   0.0677
CGKF Q > 0 4.0467   0.4115   0.5734 0.2218   0.0682   0.4862
RRR Q > 0 3.9770   0.4085   0.6456 0.2337   0.0714   0.4047
Results for the gain
Case Mean STDV
CGKF Q ≡ 0 −0.0655   0.0002   0.0027   -0.0758 0.0410   0.0286   0.1129   0.0529
RRR Q ≡ 0 0.0184   0.0077   −0.0955   −0.0038   0.2001   0.0019   0.0096   −0.0169   0.0055   0.0457 0.0020   0.0014   0.0102   0.0006   0.0216
0.0002   0.0009   0.0018   0.0006   0.0045
CGKF Q > 0 0.6249   −0.0038   −0.0017   0.4401 0.1199   0.0460   0.2033   0.1246
RRR Q > 0 0.6347   −0.0097   −0.0424   0.0111   0.0711
−0.0035   0.5409   −0.0796   0.0330   0.1417		 0.0737   0.0107   0.0487   0.0087   0.0702
0.0035   0.0886   0.1506   0.0218   0.2206

Table 1.

Comparison of simulated SMD data results of CGKF with RRR.

The gains are to be read as first column first and the second column next. For CGKF, there are only four gains associated with the two states and two measurements. But for RRR, there are ten gains associated with five states and two measurements. For the case of Q ≡ 0, all the gains should have been ideally zero but are around zero here due to the statistical percolation effect of the unforgiving noises, be it process and/or measurement. This affects not just one parameter or state but every other quantity, so the gains or any estimated quantity in the numerical algorithm can also never take their true values except perhaps with an appropriate algorithm that can capture the true values with increasing amount of data. For the Q > 0 case, the major gains marked in bold are somewhat similar.

3.3 Analysis of real airplane flight test data

This real data set discussed earlier in [6, 8, 9] is obtained along with airplane, flight test data and notations from [44]. Briefly, to explain the scenario, there is a peculiar manoeuvre when the aircraft (T 37 B) is rolling through a full rotation using the aileron, and then the elevator angle ( δ e in deg) is imparted. The coupling between the longitudinal and lateral motion is replaced by their measured values, namely the roll angle (φm), sideslip (βm), velocity (Vm), roll rate (pm), yaw rate (rm) and the angle of attack (αm). The state equations (n = 3) for the angle of attack (α), pitch rate (q) and the pitch angle (θ), respectively, are

α ̇ = − q ¯ S m V m Cos β m C L α α + C L δ e δ e + C L 0 + q + g V m Cos β m cos ( ϕ m cos ( α m cos θ + sin α m sin θ ) − tan β m p m cos α m + r m sin α m ) E28
q ̇ = q ¯   S   c ¯ I yy C m α α + C m q c ¯ 2 V q + C m α ̇ c ¯ 2 V α ̇ + C m δ e δ e + C m 0 + I zz − I xx I yy r m p m E29
θ ̇ = q cos ϕ m − r m sin ϕ m + θ 0 E30

The measurement equations (m = 4) for the angle of attack, pitch rate, pitch and normal acceleration are given by

α m = K α α − K α x α q V ; E31
q m = q ; E32
θ m = θ E33
a n m = q ¯ S mg C N α α + C N δ e δ e + C N 0 + x a n g q ̇ E34

The unknown parameters (p = 10) are C L α C L δ e C L 0 C m α C m q C m α ̇ C m δ e C m 0 θ 0 C N 0 T with the approximation C N α = C L α and C N δ e = C L δ e . The suffix δ e denotes control derivatives, and suffix zero refers to biases and all others are aerodynamic derivatives. The initial states are taken as the initial measurements and the initial parameter values are taken as (4, 0.15, 0.2, −0.5, −11.5, −5, −1.38, −0.06, −0.01, 0.2)T. At a measurement similar to the SMD system, there is an additional term Kν k which is the product of the gain matrix multiplying the innovation.

3.4 Remarks on the real flight test data results

Table 2 below compares the parameter estimates and their CRBs (in parenthesis) from the RRR [6, 8, 9], Gemson [16], (derived from the filter covariance) and CGKF (based on cost function) approaches. The parameter estimates from the first two are comparable except for the parameters C L δ e and C m q which strongly affect the airplane dynamics. However, all the parameter estimates from RRR, Gemson and CGKF are quite comparable. The CGKF estimates are within about 1σ of the RRR values as in the previous SMD case. The STDV from the CGKF (corresponding to CRB) is somewhat different from the other approaches since it follows a slightly different dynamical model than RRR or Gemson as in the SMD case.

Θ RRR Gemson CGKF
C N 0 0.2538 (0.0014) 0.2503 (0.0014) 0.2512 (0.0006)
C L 0 0.2409 (0.0021) 0.2529 (0.0018) 0.2443 (0.0019)
C L α 4.9235 (0.0164) 4.9028 (0.0168) 4.8035 (0.2199)
C L δ e 0.1554 (0.0271) 0.0879 (0.0267) 0.1653 (0.0993)
C m 0 −0.0425 (0.0009) −0.0507 (0.0024) −0.0459 (0.0001)
C m α −0.5293 (0.0079) −0.6174 (0.0211) −0.4986 (0.0345)
C m q −11.8596 (0.2402) −18.8339 (.8379) −9.3528 (2.1234)
C m α ̇ −6.8959 (0.4891) −7.1290 (1.544) −6.6730 (4.6084)
C m δ e −0.9731 (0.0177) −1.1841 (0.471) −1.0063 (0.0019)
θ 0 0.0003 (0.0021) −0.0037 (0.001) 0.0020 (0.0003)

Table 2.

Comparison of real flight test data results (Θ, σ(Θ)).

Advertisement

4. Introduction to flight data analysis of a ballistic rocket (BR)

During the development of any BR, it is necessary to carry out many flight trials and compare the flight performance with that based on pre-flight estimates. For a BR, an accurate estimation of drag coefficient is very important due to its direct impact on the system performance as it plays a very critical role in generating the firing tables. One also uses wind tunnel tests or computational fluid dynamic codes to obtain the aerodynamic characteristics. But there exist generally unavoidable errors due to wind tunnel wall interference and the limitation of wind tunnel Reynolds number. Hence, the assessment of the aerodynamic coefficient from the full-scale flight test of vehicles is an important area of activity and research. Such an analysis would help the BR as follows.

  1. If it fails en route, a real-time state estimation helps to obtain the expected impact location from range safety viewpoint.

  2. A compatibility check of measured data reduces bias and scale factor effects in the measurements. The measurement noise covariances given in the manufacturer’s catalogues being notional, such values can also be estimated.

  3. In its external ballistics, the variation of aerodynamic drag coefficient with respect to Mach number is very important.

  4. Comparison of pre-flight with flight test estimated drag coefficient helps to improve and modify the former.

4.1 State estimation of a ballistic rocket (BR)

It is possible to formulate the Kalman Filter (KF) to simultaneously estimate both the state and parameter or carry out the same sequentially. In the state estimation step, the bias, scale and the random errors are estimated, called compatibility check, and thus relatively clean data are available for parameter estimation. The BRs are generally tracked by ground based radar, which provides range, azimuth and elevation measurements. Sarkar explains in [30] the extended Kalman filter [24] together with a smoother for handling the effect of both the process and measurement noise contained in the measured flight test data. Later, it is used to estimate the aerodynamic drag coefficient (which is a parameter estimation problem) using the MMLE-OEM approach. We discuss here only his trajectory estimation by using CGKF and the reader can refer to [30] for drag estimation.

In order to track the trajectory of a BR, one can use either the dynamical or the kinematical equations. The former needs many inputs such as the forces and moments, propulsion and the control which may not be available and more so if the BR belongs to an adversary. Broadly, the three approaches all utilising the kinematic state equations for trajectory estimation with increasing accuracy are

  1. The ‘generic’ ones called α, αβ or αβγ types of filters as found in Blair [45] and Bar-Shalom and Li [46].

  2. The ‘similar’ ones like the CGKF approach which can handle similar situations.

  3. The ‘specific’ one for a given scenario like the adaptive extended Kalman filter (AEKF) such as by Gemson [16] or the ones like the adaptive limited memory filter (ALMF) as in [15].

The AEKF/ALMF deals with a specific scenario and adaptively obtains Q and R by minimising the cost function J in Eq. (16) and the steady-state gain K follows. The second CGKF handles the same specific scenario by minimising the cost function J in Eq. (17) and obtains the gain K directly. Due to the transformed unknown variables, the results for K will be somewhat different but close to AEKF. The gain being more robust, CGKF can handle similar situations. In order to account for model deficiencies or uncertainties in real cases, these constant gains can be increased from the ones based on simulated studies. The αβγ types of filters define a manoeuvre index called λ (based on a subjective choice of Q , and R, and the time between measurements) which leads to the various gains. Thus, λ being chosen subjectively, the model accountability is generic. Such filters also do not consider any cost function, whereas the second and third are two routes to tune the Kalman filter by minimising J , the normalised innovation sequence (NIS) cost function. The only way to improve the performance of αβγ types of filters is to tune the λ manually as is shown later.

The filter world kinematic model equations could consist of displacement, velocity, acceleration, jerk, slack and so on driven by inputs at the highest state derivative variables as chosen by the analyst. The inputs could be random white noise or even correlated noise. If the input is a random white noise, then the corresponding state variable and lower ones become Gauss-Markov (GM) processes of increasingly higher order.

One cannot drive the displacement by white noise since no real-world system can be instantaneously displaced due to its finite mass and moment of inertia. It is best to introduce the input at higher levels as a white or correlated Gaussian noise as in Mehrotra and Mahapatra [47], or Singer [48]. Usually, the input being a white noise acceleration, its integration provides velocity and the further integration leads to displacement. Hence, the displacement would become a third-order GM process. If the input is at the velocity level, then the displacement becomes second-order Gauss-Markov process. The input at the jerk or the slack would increase the order of the filter equations. Hence, the analyst has to choose the input acceleration at a level to provide a reasonable balance between the model order and the anticipated dynamic rate of change of the object.

4.2 Comparison of second and third order Gauss Markov system model in CGKF

Here, we consider the nine state variables as ( R , R ̇ , R ¨ , A , A ̇ , A ¨ , E , E ̇ , E ¨ ) and the three measurements as ( R , A , E ), both in polar coordinates, and the specific to real-time application in the so called PPR [30] frame. The estimation error of the state variables’ position and velocity based on the second and third order model shows it is more in the former than in the latter. This is due to the simple fact that in the second order model, the accelerations are not accounted for properly during the transition from boost phase to power off coast, when there is a rapid change in BR acceleration level, which is not taken into account in the second-order model unlike in the third order model [30].

4.3 Filter tuning using CGKF and adaptive estimation of (P0, Q , R)

We consider both the above ALMF and the simpler CGKF for real-time processing to gain confidence in the results. In the former, the choice of window length is important to reach the NIS equal to the number of measurements for filter tuning. The adaptive filter tuning of the statistics P 0 , Q and R has been carried out by varying the window length L to track the NIS Cost towards 3 as shown in Figure 1. The next Figure 2 shows the time variation of Q elements with data length of the adaptive EKF after NIS cost convergence. For a given manoeuvre in space, the choice of the coordinate system and hence the components along different axes could vary. Very rapid dynamics demand higher Q to track and slower dynamics demand lower Q. This leads to different overall constant Qs being injected in different state variables and thus the Kalman gains. In the same frame, if the origin is changed, trajectory can be hard or soft. For example, if initially the manoeuvre is very rapid in azimuth and elevation channels (with injected constant Q ), the filter cannot track the BR closely, thus giving rise to oscillatory tracking error. Other axes systems and sensitivity studies for filter statistics are available in Sarkar [30].

Figure 1.

Variation of the NIS cost with length of limited memory window.

Figure 2.

Time variation of Q elements with data length of the adaptive EKF after NIS cost convergence.

4.4 Real-time tracking using CGKF

The small differences in the gains from AEKF and CGKF are due to the duration of the transient and the steady state. At first, a set of R, A and E measurements are generated for one launch angle of the BR (45°). This data set is processed using AEKF to estimate P0, Q and R adaptively by equating the NIS cost function to 3, the number of measurement channels. After some initial transients, the Q and K elements become constant as can be seen in Figures 2 and 3, respectively. The steady state gains from the AEKF for 45° are used for processing the real data for a different launch angle of 75° of the same BR and the filter performs well. This is because the Q values from AEKF for 45 and 70°, being only slightly different, do not affect the gains in AEKF, so also in CGKF and thereby the filter performance. The NIS cost function for AEKF based on L = 5 is 3.05. For αβγ filter, the combination of λ equal to (0.002, 0.001, 0.001); (0.02, 0.01, 0.01); (0.05, 0.02, 0.02); (0.07, 0.05, 0.05) and (0.1, 0.1, 0.1) gave costs of 56.0, 17.0, 5.71, 3.03 and 2.85, respectively. These indicate the αβγ filter and the constant gain AEKF performance are close when λ = (0.07, 0.05, 0.05). Thus, the choice of gain elements from AEKF and CGKF is better than in the αβγ filter and the latter is simpler to implement.

Figure 3.

Time variation of K elements with data length of the adaptive EKF after NIS cost convergence.

Advertisement

5. Space debris re-entry

An accurate prediction of re-entry time of large orbital space debris is useful to plan hazard assessment and mitigation strategies. The database for such an analysis of large objects is the set of two line elements (TLEs) provided by agencies like USSPACECOM. The TLE sets [49] provide information regarding orbital parameters together with rate of mean motion decay and a reference parameter B* related to the ballistic coefficient B as

B ∗ = ρ o / 2 B   =   ρ o / 2 C D A eff / m E35

B represents the sensitivity of an object to air drag and B* is an adjusted value of B using the reference value of atmospheric density ρo at a reference altitude 120 km above earth. CD is the non dimensional drag coefficient, m is the mass and Aeff is the effective area of cross-section of the object. Larger B means its orbit decays faster.

5.1 Re-entry case study of US Sat. No. 25947, Soyuz

The Satellite No. 25947 is a rocket body that has been the test case for the third IADC Re-entry campaign. The sets of 72 orbital elements were made available for re-entry prediction during February 2, 2000, to March 3, 2000.

5.2 Filter-world scenario: state equations

The measurements are available in terms of the orbital parameters the semimajor axis ‘a’ and the eccentricity ‘e’ in both the simulated cases and the tracked TLE elements. The state equations governing the state variables (a, e, B) are [29]

a k = a k − 1 + Δ a k − 1 + w 1 = a k − 1 + ϕ 1 a k − 1 e k − 1 B k − 1 + w 1 E36
e k = e k − 1 + Δ e k − 1 + w 2 = e k − 1 + ϕ 2 a k − 1 e k − 1 B k − 1 + w 2 E37
B k − 1 = B k − 1 + w 3 E38

where ϕ1 and ϕ2 are the functional forms of King-Hele [50] which depend on ballistic coefficient B, ‘a’ and ‘e’, and w 1 , w 2 and w 3 are, respectively, the process noises. The subscript argument inside the brackets (∙) denotes the time instant.

5.3 Filter-world scenario: measurement equations

But, in the filter implementation process, the transformed variables, namely the predicted apogee and perigee heights, are

ha − k = a k 1 + e k and hp − k = a k 1 − e k E39

The measurements at time t(k) are of the form

ha − k = ha − k + v 1 and hp − k = hp − k + v 2 E40

where v 1 and v 2 are the measurement noises assumed to be white Gaussian with zero mean and covariance R assumed as constant. The predicted values of these heights in the state equations are updated by utilising the measured values ha and hp , respectively, the apogee height and the perigee height.

ha + hp + B + k = ha − hp − B − k + k 11 k 12 k 21 k 22 k 31 k 32 ∙ ha − ha − hp − hp − k E41

The superscripts (+) and (−) correspond to the predicted and updated values, and suffix k denotes the time instant. Further details are available in Anilkumar [26] and Anilkumar et al. [29].

5.4 Uncertainties in the state and measurement equations requiring Q and R

In general, the physical parameters like mass, shape and dimensions of the re-entry objects that vary are not available accurately. Also, the atmosphere varies randomly. Further, the tumbling effect of the body and the molecule gas surface interaction leads to uncertain and varying aerodynamic drag coefficient, which makes the prediction of re-entry time a difficult problem. The re-entry objects are mainly affected by the atmospheric drag, earth’s oblateness, solar activity index F10.7 and magnetic index Ap. However, the orbital propagator utilised in this study is a very simple model of King-Hele [50] which accounts for only the atmospheric drag effect. The present propagator assumes a mean atmospheric condition as provided by the US Standard Atmosphere [51]. This model estimates only the semimajor axis and eccentricity decay with respect to one revolution, assuming a constant scale height during one revolution. This model is sufficient for the re-entry prediction as the decay of the object is mainly governed by the air drag only. The effects of other orbital perturbations and variations in the atmospheric density are accountable through the process noise and the Kalman filter is thus able to handle it through the proper gains as will be demonstrated subsequently. In all the prediction exercises, when the semimajor axis of the object reaches a height of 120 km above the earth, it is considered to have re-entered the atmosphere. This assumption is appropriate as a reference condition since there are significant variations in the atmospheric properties above 120 km with solar, magnetic activity and local time than below this height. Also, effectively, a diffusive equilibrium predominates beyond 120 km as given in Whitten et al. [52].

5.5 Adaptive filtering approach for re-entry prediction

It was found to be adequate to obtain P0 based on the difference between the assumed initial conditions of the states (a, e) with those from the first TLE set. For state B, the deviation of initially assumed B0 with that of the derived value of B from B∗ is used. For Q and R, the heuristic estimators of Myers and Tapley [15] have been used. A careful study of data of varying length based on adaptive filtering (both by simulation and actual data) helped to assess how the estimated B, Q and R vary with data length. Recently, the RRR [6, 7, 8, 9] has found a near optimal solution for tuning the filter statistics and thus an improvement over earlier adaptive procedures.

5.6 The CGKF approach for re-entry prediction

The present study utilised the genetic algorithm (GA) in the CGKF to minimise the cost function J. The fundamentals of GA, its features and other implementation aspects can be found in [39, 40]. The values of the parameters arrived at after some trials for the present GA re-entry problem are: population size = 100; bit length = 20; probability of cross over = 0.90; probability of mutation = 0.05; number of generations for convergence = 50 and tolerance for convergence: change in cost function J between generations = 0.0001.

Starting from the 22nd TLE set, the present constant gain Kalman filter algorithm utilises a total of six gains corresponding to the three states, namely, apogee and perigee heights and the ballistic coefficient and two apogee and perigee height measurements. An important parameter in this implementation is the initial assumed value of ballistic coefficient B0. This is to be expected as the body may be tumbling, with irregular shape and with varying gas molecule and surface interaction reflected in the predicted ballistic coefficient. Further, for the drag, a very simple mean atmospheric condition is used. Figure 4 shows that as time passes, with more and more TLE data sets being available, for various initial B0 values, the predicted re-entry date comes closer. But the point is what is the best choice for the initial B0 that provides minimum variation in the predicted re-entry time right from the beginning up to the actual re-entry? This turns out to be B0 = 0.40 as shown in Figure 5. The overall problem is to find out the best possible B0 and the constant Kalman gain that predicts the re-entry time with least variation from the beginning to the end.

Figure 4.

Variation of B with time during re-entry.

Figure 5.

Variation of re-entry time with B0.

A combination of adaptive filtering and the constant gain approaches provided a set of constant gains as [0.6, 0.2, 0.2, 0.6, 0.00014 and 0.0001], as nearly optimal [26, 29]. One curious fact that may be noticed is the choice of the optimum Kalman gains. The optimal gain values for the states are larger and for B it is very small, the reason being the noise-to-signal ratio is very small for the states. Hence, the filter can track the state with a large gain value close to unity or even a small non-zero value (but not zero!). However, for the ballistic coefficient B, the gains have to be smaller in order to slowly learn from the measurements; and if these gains are larger, then the estimated B will show lot of fluctuations. The actual re-entry occurred on March 4, 2000, at 5 h 50 min. The CGKF formalism based on a mean atmosphere and approximate drag effects predicted the re-entry on March 4, 2000, at 5 h 35 min. Even the MSIS-86 model [53] could have been used. This shows once again the robustness of the constant gains has the ability to handle the inaccuracies in modelling B, as well as both the unmodelled and unmodellable state and measurement noise characteristics.

Advertisement

6. Evolution and expansion of the space debris scenario

6.1 Introduction

The evolution of the space debris scenario consisting of characterising each and every fragment can be very unwieldy. The purpose of the present study is to demonstrate that it is not necessary to follow each and every fragment in a complex environment, which demands enormous amount of computing time. It suffices to group the fragments called equivalent fragments (EQF) in the ‘a’, ‘e’ coarse bins and propagate these with time. However, the orbital and ballistic coefficients of the EQFs need to be redefined for the above purpose of time propagation in terms of the individual fragment characteristics constituting it. After time propagation, the number of fragments and their ballistic coefficient constituting the EQFs are updated based on just the measured number of individual fragments as will be explained later. This process is continued with subsequent measurements.

For studying the long-term evolution of the space debris, an initial model like in Johnson and McKnight [54], ASSEMBLE model of Anilkumar et al. [27], and Rossi et al. [55] can be assumed. At large times, the prediction could depart greatly from the real scenario due to the sensitivity of the evolution to the inaccuracies in the model parameters and the environment. There are large differences in the estimated characteristics among many debris models [26, 54]. The only way the prediction can be made to follow more closely the real situation is to update the characteristics by assimilating properly the subsequent measurements of the number density in (a, e) bins repeatedly for further evolution in time. The present procedure in addition also expands the scenario for the distribution of the ballistic coefficient of the debris as well(!) which is not generally available or measurable.

6.2 The present approach and the stochastic analog tool of Rossi et al.

The stochastic analog tool (STAT) of Rossi et al. [55] simulates the time evolution of the debris by a threefold subdivision of namely: (i) the semimajor axis ‘a’ from 6378 to 46,378 km, (ii) the eccentricity ‘e’ from 0 to 1 and (iii) the mass ‘m’ from 1 mg to 10,000 kg. The present approach considers a, log(e) and log(B) as against a, e and m of STAT. The third parameter B has been presently used because the orbital parameters are sensitive to the air drag and thus change with time. There are errors due to discretisation and approximation in specifying the arithmetic mean values for ‘a’ and geometric values ‘e’ of the EQF in the various bins. Further, there are unaccounted or even unmodellable forces during propagation. However, all such errors can be accounted for by process noise in the state equations describing the propagation of the EQF. Since the individual representative objects of each bin are propagated, the computing time is almost independent of the debris population size in both the present and STAT approaches. It is the second step that is fundamentally new and different in the present approach namely at an update apart from assimilating the measurement information it also expands the scenario to update the equivalent ballistic coefficient (EQB) for the EQFs in various bins with time.

6.2.1 Characteristics of the equivalent fragment (EQF) in terms of the fragments in a bin

Presently, with 10 divisions for each of the parameters a, e and B, a total of 1000 bins are formed. Instead of handling each and every fragment, the fragments in every bin are handled as a fewer number of EQFs. Next, to follow the dynamics of these EQFs, it is necessary to assign suitable orbital and ballistic coefficient values for these EQFs in terms of the individual fragment properties. Presently, these are set as the arithmetic mean for ‘a’, geometric mean for ‘e’ and the geometric mean also for ‘B’ of the number of fragments in each bin. As the EQFs meander across the various ‘B’ bins, their ballistic coefficients are updated. This is somewhat similar to a debris with a certain value of B moving in the atmosphere though it could change its value. A priori, how well a mean defined as above can follow the dynamics and subsequently get updated in the filter is not conceptually clear. Its adequacy can only be demonstrated from subsequent results.

6.2.2 Evolution of individual fragments as well as EQF in the bins

Initially, about 10,000 simulated debris fragments due to an explosion are considered and the later fragments due to further breakups are accounted for as source terms. The state propagation equations for both the fragments and the EQF are identical to the earlier Eqs. (38), (39) and (40). The EQF propagated based on its assigned value of suitable ‘a’ and ‘e’ and could in general end up in just within another bin. In order to redistribute the fraction of the EQFs among the bins, a heuristic rule is used as in Rossi et al. [55] that takes the ratio between the area covered by the propagated rectangle and the area of the initial rectangle as shown in Figure 6.

Figure 6.

Representation of the orbital propagation and redistribution in the eccentricity ‘e’ versus semimajor axis ‘a’ space.

Subsequently, by using the measurements of the number of debris in the bins at various times, the EQB of each EQF is updated based on the weighted average of the predicted and measured number density of the fragments in the bins. This weightage is the Kalman gain as we will see later on. Without the update of the EQB of these EQF, the filter is unable to follow the true time variation of the number density of the debris fragments in the bins. Hence, the ballistic coefficient of the EQF is aptly called as the EQB.

6.2.3 Update of the EQF characteristics

The evolution of the EQF takes place in two steps, namely: (i) the propagation of the EQFs representing all the fragments in the various bins, then, redistribution of the fragments around the adjacent bins as mentioned earlier; further breakups are also accounted for by the changed number density in the various bins, and (ii) using appropriate constant Kalman gains for obtaining an updated estimate of the number density of the fragments in the various bins and the EQB of the EQFs.

There is one subtle point in the estimation of the EQB of the EQF corresponding to various (a, e) bins. After update, the value of B for the EQB of EQF at times can fall outside the fixed bin interval. Presently, we have taken the propagation of the EQF always from the initial (a, e) condition based on the arithmetic and geometric mean, respectively. But, for a group of fragments, the above initial condition may not be the most appropriate. The EQF could have started its trajectory from anywhere inside or on the boundary of (a, e) bin whence the redistribution could have been different and thus the updated ballistic coefficient B as well. Such features arise due to the definition of the EQF characteristics and the coarseness of the bins, but one has to see if the final results are meaningful and acceptable.

6.2.4 Real-world (individual fragments) and filter-world (EQF) scenarios

The state and measurement equations in the real-world and filter-world scenarios are given in Table 3. In the filter state equations, the binning, formation of EQF, propagation and redistribution all lead to modelling error and need process noise to handle the situation. In a real-world scenario, there would be measurement noise due to inaccuracies in the assigned orbital characteristics of the individual fragments. In simulation studies, the propagation of each and every one of the individual debris fragments and counting their number in the various bins lead to no measurement noise.

Quantity Real World Scenario for each fragment Filter World Scenario for each EQF
The state variables The (a, e, B) of each fragment. The (a, e, B) of each EQF.
The initial conditions The initial (a, e) of each fragment. For the EQF from anywhere in the (a, e) bin.
The state input Complex environment. Only the air drag effect is considered.
The state process Random variations in the real environment. The inaccuracies in assigning (a, e, B), binning, its propagation, redistribution and the environment.
The measured variables The number of individual fragments in the various bins. The EQFs are propagated with only air drag and later converted to the number of objects in each of the bins.
The measurement noise Measurement errors due to tracking and data processing. No measurement noise as the EQFs are propagated and using the changed values of (a, e) are assigned to appropriate bins.

Table 3.

Real world without binning and filter world with binning (for simulation studies).

The uncertainties in the initial EQB values of the EQF in the state equations (shown in Table 3) are improved by the filter by using a certain length of data. In the present study, the constant Kalman gains have been derived as explained later.

6.3 The present constant gain Kalman filter approach

From simulated studies, the number of debris fragments in each three-dimensional (a, e, B) bin is known exactly. The Kalman filter by using the constant gains and the updated number of objects at various times is able to track closely the true number of fragments. Similarly, the measurements can be assimilated and the scenario expanded to get the EQB.

6.3.1 Filter-world state equations

Thus, the states presently considered in every one of the (a, e) bins are the number of objects N and their EQB. The (a, e) bins are not changing and the EQF moves in the (a, e) plane like any other single fragment and later gets redistributed based on a certain rule. The various EQFs are specified by: (i) their number in each (a, e) bin; (ii) the equivalent semimajor axis, (iii) the equivalent eccentricity and (iv) the EQB.

The state equations for the EQF in the various bins between measurements are

d N / dt = Σ propagation across the a e bins and redistribution + source terms + state noise E42
d B / dt = 0 + state noise = 0 E43

6.3.2 Filter-world measurement update equations

The true number density NM in the various bins is obtained in simulation by propagating each and every individual debris fragment. This is used to update the predicted number density based on propagated and redistributed EQF as well as update the ballistic coefficient of the EQF in the various bins as given by

N + = N − + K N . N M – N − ; E44
B + = B − + K B . N M – N − E45

with the pre- and post-updated values denoted, respectively, by the superscripts (−) and (+). The gains KN and KB, respectively, correspond to the number density and the equivalent ballistic coefficient.

Presently, the number of constant gains KN and KB to be estimated is 200 with two for each of the 100 (a, e) bins. These are obtained based on minimising the cost function.

J = 1 / N Σ ν k R − 1 ν k T E46

Σ denotes the summation over all bins and times. The constant Kalman gains were obtained by using the genetic algorithm [39, 40]. The different parameters used in GA implementation are as follows: population size = 200; bit length = 20; probability of cross over = 0.90; probability of mutation = 0.05; Convergence: number of generations 50 or alternately change in J between generations 0.0001.

The whole of Kalman filter process can be summed up in a simple way. One can have the evolution of the state (without knowing how) generated by any random process. The time variation of the state can even be assumed to be given. The measurements could be noisy or even exact. In order to track the state and also follow it smoothly by reducing the fluctuations, a simple filter can be designed with the states remaining constant between measurements. For zero Kalman gain, the filter will learn nothing from the measurements and the state will remain at the initial values. For unity Kalman gain, the state will follow the measurements. In between, there is a range of gain for which the difference between the predicted state and the measurement is minimised in a suitable sense over the range of data. For slow and fast state dynamics, gains near zero and unity, respectively, would be appropriate. Further, if the random process is known to have an inaccurate or unknown parameter, they can also be handled by additional constant gains.

6.3.3 Evolution of debris objects generated due to explosions

A single explosion at a typical altitude of 800 km and eccentricity 0.00045 resulting in about 10,000 debris of varying ballistic coefficients is simulated using the ASSEMBLE model [26, 27]. These objects were propagated accounting for only the atmospheric drag effect for a period of 600 days and thus generate the (a, e, B) data of the objects. Among many studies, updating both the number of objects and the EQB gave the best results and is described here for brevity. Further experiments like initial explosion followed by one breakup, two break ups and some launch activities were carried out. The filtering process reduces the errors in the estimates, but not below a certain value due to continuous occurrence of error due to binning, propagation and redistribution leading to a non-zero K. In all the subsequent figures, the symbol (o) denotes the true, (solid line) filter and (dashed line) with no update using K = 0.

6.3.4 Evolution of a single breakup and additional debris

Figure 7 provides the results by updating both the number density and the EQB for the typical B bin (0.6056, 1.6476). Figure 8 shows the variation of the constant Kalman gain KN across the semimajor axes bins for four ballistic bins are always between zero and unity since the state, namely the number density, is measured.

Figure 7.

Breakup evolution in B bin (0.6056, 1.6476).

Figure 8.

Kalman gains for number of objects (KN).

Figure 9 shows the KB for various values of the semimajor axes bins. However, this takes positive and negative values. Such a thing can happen since the ballistic coefficient though a state has not been measured. Further, in other experiments that were performed, it was noted that the estimated ballistic coefficient with time in typical bins is generally within the limits of the ballistic coefficient bin values. However, at times, they move somewhat outside the limits of the bin values. The initial condition for EQF propagation from the mean values of the bins, though it could have started from anywhere inside the bins, the subsequent propagation and redistribution error could be responsible for such a behaviour. It is best to accept the approach here as the ability to mimic the dynamical behaviour.

Figure 9.

Kalman gains for ballistic coefficients (KB).

At the beginning 10,000 fragments were introduced and subsequently an additional 300 fragments were introduced after 120 days very much like the real-world scenario where the debris are growing but not too rapidly. Even after adding the new source terms, the constant Kalman gains obtained based on the initial cloud evolution have been used for further evolution. In general, the constant gains are robust around a range of the estimated values. Hence, even if the subsequent results are non-optimal, they are adequate to obtain acceptable estimates.

6.3.5 Evolution of a single breakup followed by more than one breakup and some launch activities

To simulate the real-world scenario, two explosions and some launches are introduced during the evolution. Once again, the constant gains obtained using the primary debris clouds suffice for all later cases as well! The results provided in Figure 10 show that the present model is able to track the number of objects even in such evolution process.

Figure 10.

Estimated and observed number of objects. (Two explosions and some new launches).

6.3.6 Application to a typical real-world scenario

The catalogued TLE data of 335 debris objects in near circular orbits in the perigee and apogee bin of 700–800 km from October 1998 to September 1999 were chosen, assuming an initial ballistic coefficient for the EQFs is propagated and updated using first eight observations. A constant eccentricity for all the EQFs in all the bins is assumed since the bin size is just 10 km (unlike in the simulated scenario where it was 150) km but their semimajor axis corresponds to the mid-value of the bin. The Kalman gains from simulation studies were once again used to analyse the real-world data(!), thus demonstrating the robustness of the constant Kalman gains. The innovation here is given by the difference between the TLE data and the predicted number density. The 335 objects observed in the apogee/perigee bin from October 1998 were tracked for the next 12 months to obtain the number of objects in the semimajor axis bins. By tracking the same objects, Figure 11 provides their number for the 12 months in the 10 different semimajor axis bins.

Figure 11.

Variation of the number density of the debris in the bins with time.

Figure 12 shows a comparison of the number of objects from the present approach with that observed for the semimajor axis bin of (7150, 7160) km. Considering 10 equivalent objects rather than propagating and monitoring all of the 335 objects, the match is quite good.

Figure 12.

Comparison of the estimated time-varying debris number density with TLE values.

Advertisement

7. Conclusions

The present CGKF approach has been demonstrated with many examples and in particular the evolution of thousands of space debris fragments. This formalism can be used even in massive atmospheric data assimilation and weather prediction problems that have tens of thousands of states and measurements.

Advertisement

Acknowledgments

The author sincerely thanks all collaborators referred to in this chapter. It has been very kind of Prof. Felix Govaers to have invited me to contribute a chapter in this book and also to INTECHOPEN for publishing it.

References

  1. 1. Kalman RE. A new approach to linear filtering and prediction problems. Transactions of the ASME - Journal of Basic Engineering. 1960;82(Series D):35-45
  2. 2. Brown R, Hwang P. Introduction to Random Signals and Applied Kalman Filtering, with MATLAB Exercises. 4th ed. John Wiley and Sons; 2012
  3. 3. Kailath T. An innovations approach to least-squares estimation part I: Linear filtering in additive white noise. IEEE Transactions on Automatic Control. 2016;13(6):646-655
  4. 4. Klein V, Morelli EA. Aircraft System Identification: Theory and Practice. AIAA Edn. Series; 2006
  5. 5. Jategaonkar RV. Flight Vehicle System Identification: A Time Domain Methodology. Vol. 216. AIAA; 2006
  6. 6. Shyam MM, Naik N, Gemson RMO, Ananthasayanam MR. Introduction to the Kalman filter and tuning its statistics for near optimal estimates and Cramer Rao bound, TR/EE2015/401; IIT Kanpur; 2015. Available from: http://arxiv.org/abs/1503.04313
  7. 7. Ananthasayanam MR, Shyam Mohan M, Naik N, Gemson RMO. A heuristic reference recursive recipe for adaptively tuning the Kalman filter statistics part-1: Formulation and simulation studies. Sadhana. 2016;41(12):1473-1490
  8. 8. Shyam Mohan M, Naik N, Gemson RMO, Ananthasayanam MR. A heuristic reference recursive recipe for adaptively tuning the Kalman filter statistics part-2: Real data studies. Sadhana. 2016;41(12):1491-1507
  9. 9. Ananthasayanam MR. A reference recursive recipe for tuning the statistics of the Kalman filter. In: de Oliveira Serra GL, editor. Kalman Filters-Theory for Advanced Applications, INTECHOPEN. DOI: 10.5772/intechopen.71961
  10. 10. Narasimha R. Performance reliability of high maintenance systems. Journal of the Franklin Institute. 1975;303:15-29
  11. 11. Ananthasayanam MR. A relook at the concepts and competence of the Kalman filter. In: Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit; 5-8 January 2004; Reno, Nevada. AIAA-2004-571
  12. 12. Ananthasayanam MR. Fascinating perspectives of state and parameter estimation techniques, AIAA-2000-4319
  13. 13. Sorenson HW. Least squares estimation from Gauss to Kalman. IEEE Spectrum. 1970;7:63-68
  14. 14. Mehra R. Approaches to adaptive filtering. IEEE Transactions on Automatic Control. 1972;17:903-908
  15. 15. Myers KA, Tapley BD. Adaptive sequential estimation with unknown noise statistics. IEEE Transactions on Automatic Control. 1976;21:520-525
  16. 16. Gemson RMO. Estimation of aircraft aerodynamic derivatives accounting for measurement and process noise by EKF through adaptive filter tuning [PhD thesis]. Bangalore, India: Department of Aerospace Engineering, Indian Institute of Science; 1991
  17. 17. Gemson RMO, Ananthasayanam MR. Importance of initial state covariance matrix for the parameter estimation using adaptive extended Kalman filter. In: Proceedings of AIAA Atmospheric Flight Mechanics, AIAA-98-4153. 1998. pp. 94-104
  18. 18. Wilson KC. An optimal control approach to designing constant gain filters. IEEE Transactions on Aerospace and Electronic Systems. 1972;8(6):836-842
  19. 19. Cook G, Dawson DE. Optimum constant-gain filters. IEEE Transactions on Industrial Electronics and Control Instrumentation, IEEE. 1974
  20. 20. Grimble MJ, Jukes KA, Goodall DP. Nonlinear filters and operators and the constant-gain extended Kalman filter. IMA Journal of Mathematical Control and Information. 1984;1:359-386
  21. 21. Kobayashi T. Application of a constant gain extended Kalman filter for in-flight estimation of aircraft engine performance parameters, NASA/TM—2005-213865. 2005
  22. 22. Liu X, Yu Y, Li Z, Iu HHC, Fernando T. A novel constant gain Kalman filter design for nonlinear systems. Signal Processing. DOI: 10.1016/j.sigpro.2017.01.002
  23. 23. Song TL, Speyer JL. A stochastic analysis of a modified gain extended Kalman filter with applications to estimation with bearing only measurements. IEEE Transactions on Automatic Control. 1985;(10):940-949
  24. 24. Gelb A. Applied Optimal Estimation. Cambridge, Massachusetts, MIT Press; 1974
  25. 25. Sugiyama N. System identification of jet engines. Journal of Engineering for Gas Turbines and Power. 2000;122:19-26
  26. 26. Anilkumar AK. New perspectives for analyzing the breakup, environment, evolution, collision risk and reentry of space debris objects [PhD thesis]. India: Department of Aerospace Engineering, Indian Institute of Science; 2004
  27. 27. Anilkumar AK, Ananthasayanam MR, Subba Rao PV. A posterior semi-stochastic low earth debris on-orbit breakup simulation model. Acta Astronautica. 2005;57:733-746
  28. 28. Ananthasayanam MR, Anilkumar AK, Subba Rao PV. New approach for the evolution and expansion of space debris scenario. Journal of Spacecraft and Rockets. 2006;43(6):1271-1282
  29. 29. Anilkumar AK, Ananthasayanam MR, Subba Rao PV. A constant gain Kalman filter approach for the prediction of re-entry of risk objects. Acta Astronautica. 2007;61:831-839
  30. 30. Sarkar AK. New perspectives in state and parameter estimation of flight vehicles for off line and real time applications [PhD thesis]. Bangalore, India: Department of Aerospace Engineering, Indian Institute of Science; 2004
  31. 31. Viswanath. Study of constant gain Kalman filter for basic systems with preliminary study of ISAR problem [MTech thesis]. India: IIT Kanpur; 2018
  32. 32. Philip NK. A study of the estimation schemes with controller for the final phase of a rendezvous and docking mission of spacecraft [PhD thesis]. Bangalore, India: Department of Aerospace Engineering, Indian Institute of Science; 2001
  33. 33. Philip NK, Ananthasayanam MR. Relative position and attitude estimation and control schemes for the final phase of an autonomous docking mission of spacecraft. Acta Astronautica. 2003;52:511-522
  34. 34. Anand Raj R. Kalman filter estimation of ionospheric tec and differential instrumental biases over low latitude using dual frequency GPS observations [MSc thesis]. Bangalore, India: Department of Aerospace Engineering, Indian Institute of Science; 2006
  35. 35. Anandraj R, Ananthasayanam MR, Mahapatra PR, Soma P. Estimation of Low Latitude Ionospheric TEC from Dual Frequency GPS Observations by Using Adaptive Kalman Filtering Technique. Chicago: International Union of Radio Science. Union Radio Scientifique Internationale. XXIX General Assembly; 2008
  36. 36. Basil H. Studies on a low cost integrated navigation system using MEMS-INS and GPS with adaptive and constant gain Kalman filters [MSc thesis]. Bangalore, India: Department of Aerospace Engineering, Indian Institute of Science; 2005
  37. 37. Basil H, Ananthasayanam MR, Puri SN. An application of constant Kalman gain approach for the problem of integrating the MEMS-INS and the GPS, JS-65. In: KSAS-JSASS Joint Symposium; Seoul, Korea; 2004
  38. 38. Yadav A, Naik N, Ananthasayanam MR, Gaur A, Singh YN. A constant gain Kalman filter approach to target tracking in wireless sensor networks. In: Industrial and Information Systems (ICIIS), 2012 7th IEEE International Conference on. pp. 1-7
  39. 39. Deb K. Optimization for Engineering Design: Algorithms and Examples. Prentice-Hall of India; 1995
  40. 40. Goldberg DE. Genetic Algorithms in Search, Optimization & Machine Learning. Addison-Wesley; 2000
  41. 41. Billard R. Aryabhata and Indian Astronomy. Indian Journal of History of Science. 1977;12(2):207-224
  42. 42. Billard R. L’astronomie Indienne. Paris: Publications De Ecole Francaise D’Extreme-Orient; 1971
  43. 43. Sarma KV. Tradition of Aryabhatiya in Kerala: Revision of planetary parameters. Indian Journal of History of Science. 1977;12(2):194-199
  44. 44. Maine RE, Iliff KW. Programmer’s manual for MMLE3, a general Fortran program for maximum likelihood parameter estimation, NASA TP-1690. 1981
  45. 45. Blair WD. Fixed-gain, two-stage estimators for tracking maneuvering targets. Report No. NSWCDD/TR-92/297. 1992
  46. 46. Bar-Shalom Y, Rong Li X, Kirubarajan T. Estimation with Applications to Tracking and Navigation, Theory, Algorithm and Software. John Wiley and Sons; 2000
  47. 47. Mehrotra K, Mahapatra PR. A jerk model for tracking highly maneuvering targets. IEEE Transactions on Aerospace and Electronic Systems. 1997;33(4):1094-1105
  48. 48. Singer RA. Estimating optimal tracking filter performance for manned maneuvering targets. IEEE Transactions on Aerospace and Electronic Systems. 1970;6:473-483
  49. 49. Web sites for two line element (TLE) sets. Available from: (http://www.geocities.com/iss25544/#TLEs) or (www.wingar.demon.co.uk) or (http://www.celestrak.com)
  50. 50. King-Hele D. Satellite Orbits in an Atmosphere: Theory and Applications. London: Blackie; 1987
  51. 51. NOAA, NASA and USAF, U.S. US Standard Atmosphere. Washington, DC: Government Printing Offices; 1976
  52. 52. Whitten RC, Vaughan WW. Guide to Reference and Standard Atmosphere Models, AIAA-G-003-1988. 1990
  53. 53. Hedin AE. MSIS-86 thermospheric model. Journal of Geophysical Research. 1987;92:4649-4662
  54. 54. Johnson NL, McKnight DS. Artificial Space Debris. Orbit Book Company; 1987
  55. 55. Rossi A, Cordelli A, Farinella P, Anselmo L. Collisional evolution of the earth’s orbital debris cloud. Journal of Geophysical Research. 1994;99(E11):23195-23210

Written By

Mudambi R. Ananthasayanam

Submitted: 05 March 2018 Reviewed: 02 October 2018 Published: 23 November 2018

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Continue reading from the same book

View All
Chapter 4 Statically Fused Converted Measurement Kalman Filt...

By Gongjian Zhou and Zhengkun Guo

1368 downloads

Chapter 5 A Scalable, FPGA-Based Implementation of the Unsce...

By Jeremy Soh and Xiaofeng Wu

1633 downloads

Chapter 6 Novel Direct and Accurate Identification of Kalman...

By Rajamani Doraiswami and Lahouari Cheded

1299 downloads