Comparison of simulated **SMD** data results of **CGKF** with **RRR.**

## 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

where ‘**x**’ and ‘Z’ are, respectively, the state and measurement equations of size (n × 1) and (m × 1); **u** is the control input 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 **EKF**) formulation can be written as

or equivalently

where ‘**u**’ is not shown for brevity. The formal solution for the above filtering problem can be summarised following Brown and Hwang [2] as

We assume that

with **f**/∂**X**) evaluated at **h**/∂**X**) is evaluated at

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

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 **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

where **P0**, **Q** and **R** varies with time. The estimation of the system parameters **RRR** [6, 7, 8, 9] or the heuristic approach of Myers and Tapley [15] for

### 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

gets simplified to determine only the constant gain matrix **K**(k) by subsuming **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 **K** can be obtained by making the above normalised innovation cost function equal to the number of measurements by assuming the above **MMLE-OEM** [4, 5]. There could be some differences in the gain values obtained from the adaptively or manually tuned **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.

## 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.

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.

## 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

where **Θ = [****]**^{T} has the true value **Θ**_{true} **= [4, 0.4, 0.6].** **X** = [**]**^{T} of size [(n + p) × 1] which in this case is (5 × 1). The measurement equation is given by.

where **H** =

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

where **K**. These along with the parameter vector **Θ** have to be estimated by minimising the earlier mentioned innovation cost function **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 [**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 |

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 (_{m}), sideslip (β_{m}), velocity (V_{m}), roll rate (p_{m}), yaw rate (r_{m}) 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

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

The unknown parameters (p = 10) are ^{T}. At a measurement similar to the **SMD** system, there is an additional term

### 3.4 Remarks on the real flight test data results

Table 2 below compares the parameter estimates and their **CRB**s (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 **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.

CGKF | |||
---|---|---|---|

0.2538 (0.0014) | 0.2503 (0.0014) | 0.2512 (0.0006) | |

0.2409 (0.0021) | 0.2529 (0.0018) | 0.2443 (0.0019) | |

4.9235 (0.0164) | 4.9028 (0.0168) | 4.8035 (0.2199) | |

0.1554 (0.0271) | 0.0879 (0.0267) | 0.1653 (0.0993) | |

−0.0425 (0.0009) | −0.0507 (0.0024) | −0.0459 (0.0001) | |

−0.5293 (0.0079) | −0.6174 (0.0211) | −0.4986 (0.0345) | |

−11.8596 (0.2402) | −18.8339 (.8379) | −9.3528 (2.1234) | |

−6.8959 (0.4891) | −7.1290 (1.544) | −6.6730 (4.6084) | |

−0.9731 (0.0177) | −1.1841 (0.471) | −1.0063 (0.0019) | |

0.0003 (0.0021) | −0.0037 (0.001) | 0.0020 (0.0003) |

## 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.

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

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.

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

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 **BR**s 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

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

The ‘similar’ ones like the

**CGKF**approach which can handle similar situations.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 **K** follows. The second **CGKF** handles the same specific scenario by minimising the cost function **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 **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 (**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 **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 **Q**s 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].

### 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.

## 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 (**TLE**s) 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** 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. **C**_{D} is the non dimensional drag coefficient, **m** is the mass and **A**_{eff} 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]

where ϕ1 and ϕ2 are the functional forms of King-Hele [50] which depend on ballistic coefficient **B**, ‘a’ and ‘e’, and

### 5.3 Filter-world scenario: measurement equations

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

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

where **R** assumed as constant. The predicted values of these heights in the state equations are updated by utilising the measured values

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 F_{10.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.

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.

## 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 **EQF**s 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 **EQF**s 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 **EQF**s 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 **EQF**s, it is necessary to assign suitable orbital and ballistic coefficient values for these **EQF**s 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 **EQF**s 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.

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 **EQF**s 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 **EQF**s.

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. |

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

#### 6.3.2 Filter-world measurement update equations

The true number density **N**_{M} 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

with the pre- and post-updated values denoted, respectively, by the superscripts (−) and (+). The gains **K**_{N} and **K**_{B,} respectively, correspond to the number density and the equivalent ballistic coefficient.

Presently, the number of constant gains **K**_{N} and **K**_{B} to be estimated is 200 with two for each of the 100 (a, e) bins. These are obtained based on minimising the cost function.

**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 **K**_{N} 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 9 shows the **K**_{B} 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.

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.

#### 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 **EQF**s 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 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.

## 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.

## 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.