1. INTRODUCTION
In a review of uninhabited surface vehicles (USVs) (also referred to as unmanned surface vehicles) by Motwani (Reference Motwani2012), it can clearly be seen that such craft are now being used in a plethora of marine related applications. For an example, in June 2011, the Ministry of Land, Transport and Maritime Affairs, South Korea, announced a four year USV development programme worth a total of $18.5million (Martin, Reference Martin2012). Missions envisaged include surveillance, research and monitoring of the oceans. In particular, these vehicles will be used for exclusive economic zone protection duties with a remit for policing illegal fishing and military border intrusions by North Korea.
Irrespective of its allotted task, however, all USVs have two common features. Firstly, such vehicles require robust navigation, guidance and control systems in order to cope with possible changes in the dynamic behaviour of a vehicle which may occur owing to the deployment of different payloads, amendments to mission requirements and varying environmental conditions; and secondly, unlike large commercial ships and warships that are equipped with high specification navigational aids such as radio beacons, radar, gyroscopic compasses and inertial measurement units (IMUs), the navigation suites for USVs are invariably low-cost.
As reported by Sharma et al. (Reference Sharma, Naeem and Sutton2012), in recent years at Plymouth University, the Springer USV has been designed, built, and continues to be developed. The work presented was concerned with the development of a sophisticated autopilot for Springer, whereas this paper explores the potential application of interval Kalman filtering (IKF) in the design of a navigation system for the vehicle. The investigation was instigated as the IKF has received minimal attention for application in such roles, even though, as will be discussed, its fundamental quality of being able to compute rigorous bounds to the optimal estimate can be advantageous when the model parameters themselves are uncertain.
The structure of this paper is as follows: section 2 details the state-space yaw model of the Springer used throughout this paper; section 3 then briefly discusses traditional Kalman filtering (KF), with particular attention to its use in surface navigation systems, and evidences its limitations when the system model is imprecise; an introduction to IKFs is given in section 4, and its application to the Springer yaw model shown in section 5, from which certain characteristics are gleaned leading to a discussion and review of techniques being developed with regard to these. Concluding remarks are in section 6.
2. THE SPRINGER YAW DYNAMICS AND NAVIGATION SUITE
The yaw dynamics of the Springer USV can be represented as a second-order state-space model by Naeem et al. (Reference Naeem, Xu, Sutton and Tiano2008) as shown below. The propulsion system is based on two trolling motors, one situated on each hull. If n 1 and n 2 represent the rpm of each motor (Figure 1), then while n 1=n 2 the vessel moves in a straight line, whereas differences in the two revolution rates enable the vessel to steer.
Defining n c and n d as
it is clear that steering is controlled by varying n d while keeping n c constant. Several trials have been carried out at Roadford Reservoir in Devon, UK, in which the vehicle was driven for some calculated manoeuvres, maintaining a constant n c of 900 rpm. During these trials, both the differential thrust applied to the motors, n d, and the heading angle of the vehicle, obtained from on board compasses, were recorded. System identification (SI) techniques were then applied and the following state-space model of the yaw dynamics was obtained and validated (Naeem et al. Reference Naeem, Xu, Sutton and Tiano2008):
with
and a sampling time of 1 s, where u(k) represents the differential thrust input in rpm and y(k) the heading angle in radians.
While full details of the Springer's hardware can be found in the Journal of Navigation (Sharma et al., Reference Sharma, Naeem and Sutton2012), for the purposes of this paper it is sufficient to mention that the sensor suite combines a Global Positioning System (GPS), three different types of digital compasses, a speed log, and a depth sensor, interfaced to a PC on board the vessel via a NI-PCI 8430/8 (RS232) serial connector.
3. A BRIEF OVERVIEW OF KALMAN FILTERING IN USV SURFACE NAVIGATION
Since its inception in 1960 and its initial application to spacecraft navigation, the KF has been used extensively in innumerable applications. While a detailed study of the KF can be found in Chui and Chen (Reference Chui and Chen2008), the main results are summarised in Appendix A.
The algorithm's inherent structure allows it to naturally combine measurements from various sensors, weighing their respective precisions. This has prompted the use of the KF as a tool for fusing data from low-cost sensors to obtain synergistically highly reliable estimates that would otherwise require more precise sensors. Whereas in spacecraft attitude estimations, KFs are developed for highly accurate inertial sensors, multi-sensor data-fusion has been prominent in the case of land and sea surface-vehicle navigation. Particularly, the integration of low-cost inertial navigation systems (INS) with GPS localisation data has been a commonly adopted strategy in these vehicle types. For example, the USV being developed at Virginia Tech University (VaCAS, 2011), uses differential GPS (DGPS) together with a puck-sized micro-electro-mechanical system (MEMS) based inertial sensor offered by MicroStrain Inc that is popularly used in mobile robotic applications (MicroStrain, 2012).
Others have used KFs to combine GPS, IMU, as well as magnetic compass sensor readings. For example, Zhang et al. (Reference Zhang, Gu, Milios and Huynh2005) describe the use of an unscented KF (UKF) to combine a low-cost IMU, GPS and digital compass using a sophisticated dynamical model of the USV. Others have successfully implemented KF-based USV navigation without IMUs altogether. In the case of Springer (section 2), data from the digital compasses are combined using various data-fusion architectures based on KFs (Xu, Reference Xu2007). The use of redundant data (by using three separate compasses simultaneously) allows for the construction of fault-tolerant navigation systems. Successful sea trials have demonstrated the effectiveness of the navigation systems of Springer. Another example is the USV Charlie, equipped solely with GPS and magnetic compass, which uses an extended KF (EKF) (Caccia et al., Reference Caccia, Bibuli, Bono and Bruzzone2008). The reader interested in the workings of the UKF and EKF may consult Aich and Madhumita (Reference Aich and Madhumita2010).
The basic KF scheme yields a statistically optimal estimate only for linear systems with white Gaussian system and measurement noises with known covariances, when an accurate description of the deterministic matrices, the estimated initial state, and its covariance are available. When these hypotheses are not met, the effectiveness of the KF can be compromised. As an example, consider a computer simulation of a KF estimate of the Springer heading for which the previously described Springer yaw dynamics model (Equations 2 to 4) is used for the predictive stage of the filter, while measurements of the heading angle are obtained directly from one of the compasses and used in the corrective stage. Although dynamic models for each of the compasses were obtained by Xu (Reference Xu2007), since these are several orders of magnitude faster than the Springer yaw dynamics, in practice they can be assumed to provide the instantaneous heading of the vehicle, to within a certain accuracy and precision. Assuming that the system state is affected by a random input disturbance that follows a zero-mean white Gaussian noise sequence with covariance matrix Q=diag{1,1}×10−9, and that compass readings provide accurate but imprecise measurements that likewise can be described as a white Gaussian noise sequence, with mean equal to the true heading, in degrees, and a standard deviation of 2°, the state and measurement models for the yaw dynamics of Springer are given by:
and
with
y(k) being the compass measurement in degrees. In order to illustrate the effect of deficient system modelling in the quality of the KF estimate, consider that the actual Springer yaw dynamics differ slightly from this model: specifically, that the actual value of the first component of the A matrix is 1% less than the modelled value, all other values being accurately modelled. Then, based on simulated values of ω(k) and υ(k) and the input values of n d shown in Figure 2a, the KF estimate obtained is shown against the true heading in Figure 2b. For the sake of accentuating the effect of the erroneous modelling in the KF estimate, a second KF estimate is shown for which the true model of the system was assumed.
The traditional KF relies on known statistical models (Q≡cov(ω), R≡cov(υ)) to describe uncertainty. However, as suggested earlier, these a priori statistics may not be known accurately, or be changing over time. For example, GPS accuracy is affected by the positions of the satellites, interference of the radio signal, physical barriers to the signal like mountains or those due to atmospheric conditions. Compass readings may be affected by interfering electromagnetic fields. In this scenario, the measurement noise covariance R must be continuously adapted to accommodate the changes in sensor accuracy in order to achieve good performance from the KF. Likewise, the varying sea roughness and other random effects that affect the system should ideally be contemplated in a varying system noise covariance Q that most accurately reflects the conditions at each moment. The KF also assumes that the deterministic dynamics described by matrices A, B and C are modelled precisely; however, changing payload, mean wind speeds or current forces, etc would translate into slowly-varying dynamics that would require the model to be continuously updated. While there are many techniques which try to increase the KF's reliability by incorporating adaptive mechanisms, such as online modification of the covariance matrices based on fuzzy logic (Abdelnour et al., Reference Abdelnour, Chand and Chiu1993; Xu et al., Reference Xu, Chudley and Sutton2006), such techniques cannot in general guarantee the virtue of the estimates in all situations, or provide rigorous error bounds.
4. INTERVAL KALMAN FILTERING
Interval analysis is a field of mathematics that began to be formally studied in the 1950s with the intention of finding a way to bound rounding errors in finite precision numerical computations. As computer representation is limited by the machine epsilon, only a small subset of real numbers can be represented accurately on a computer, but every real number can be represented by an enclosure consisting of two bounding, machine-representable values, by appropriate upward and downward rounding. Consequently, an interval arithmetic (IA) was defined (Moore, Reference Moore1966) in such a way that calculations on interval values yield further intervals which guarantee to enclose (though not necessarily be equal to) the range of possible results, a fundamental quality known as the inclusion property. Hence, if an initial bound for the error is known, then the error propagation inherent in any calculation carried out with finite precision is automatically bounded by the result of computing this propagation using IA. Since the publication of the first book on interval analysis (Moore, Reference Moore1966), there has been a keen interest in applying interval-based verified-computing in numerous fields, ranging from computer-assisted proofs in mathematical analysis (Lanford III, Reference Lanford1986) to practical engineering and industrial applications (Corliss, Reference Corliss and Ullrich1990). An introduction to interval analysis with a review of applications can be found in Kearfott (Reference Kearfott1996) and in Alefeld and Mayer (Reference Alefeld and Mayer2000).
In the late 1990s a KF applied to dynamical interval systems, that is, systems whose parameters are described in terms of intervals, was proposed by Chen et al. (Reference Chen, Wang and Shieh1997). Uncertainty in system modelling is often naturally described in this form. Physical parameters are usually not known precisely but specified with a certain tolerance, or have a varying nature: for example, as indicated earlier, the dynamics of marine surface vehicles depend, among other parameters, on their mass, which varies depending upon the number of passengers aboard. Even in the case of USVs, payload may vary depending on the mission as on board equipment may be configurable, and fuel-driven USVs typically have a large fuel capacity to total platform weight ratio, making such variations significant. In the Springer state-space model (Equations 2 to 4), the coefficients of the matrices were obtained using SI techniques using a certain data set registered during a specific trial. During other trials, the values obtained may vary slightly due to reasons such as those previously outlined. The effect of incorrectly describing these values was illustrated in Figure 2. However, if a whole range of trials is performed, the varying results can be enclosed in intervals, resulting in an interval system model that contains each of the individual models obtained. In this way, all the possible dynamics are taken into account.
Chen et al. (Reference Chen, Wang and Shieh1997) used the theory of IA to construct a KF for the interval system model. They obtained the IKF equations (Figure 3) using the same derivation as the regular KF. Suppose some elements of the matrices A, B and C are uncertain within some definite bounds. The system can then be described by:
where MI=M±∆M=[M−|∆M|, M+|∆M|] for M∈{A,B,C}, and ω(k) and υ(k) are white noise sequences with zero-mean Gaussian distributions with known covariances cov(ω)=Q, cov(υ)= R, and E[ω(l) υT(k)]=0 ∀l,k, E[x(0) ωT(k)]=0, E[x(0) υT(k)]=0∀k.
The IKF algorithm is summarised by the equations shown in Figure 3 (Chen et al., Reference Chen, Wang and Shieh1997), which mimic those of the ordinary KF but are described in terms of intervals. Given an initial estimate ${\bi \hat x}^{\bf I} (0)$ and its uncertainty, characterized by ${\bi P}^{\bf I} (0) \equiv {var} [{\bi \hat x}^{\bf I} (0)]$, together with the input to the system and the output measurement at each time-step, the resulting state estimate is an interval vector ${\bi \hat x}^{\bf I} (k)$ at each time-step k, providing an upper and lower boundary to the estimate, as illustrated in Figure 4.
Having been derived from the same principles, the IKF is statistically optimal in the same sense as the standard KF, and it maintains the same recursive formulation. However, the main advantage of the computed interval estimates, as opposed to point estimates, is that they guarantee to contain all the KF estimates of the individual models contained in the interval model, consequence of the inclusion property of IA by which if initial imprecise data is enclosed within rigorous bounds, then computation with these bounds will carry on yielding rigorous bounds of the actual solution range.
This can be important if one is to have confidence in the estimate. If the true system dynamics are known to be contained in the interval model, then the IKF provides a guaranteed enclosure of an optimal state estimate. While the precise value of this estimate will not be known, an interval may be acceptable for the required purpose, for example, if the object is to maintain a state variable between two limiting values, as long as the interval estimates remain within these limits no control action is required. Likewise, should the estimation boundaries permeate into an undesired operating region, this can be used to raise an alarm or trigger some other contingency mechanism.
5. IKF SIMULATION AND DISCUSSION
To illustrate the application of the IKF algorithm to the estimation of the yaw angle of the Springer USV, consider again the yaw dynamics modelled by Equations (5) to (7), but that the values of B are given with a tolerance of 25%. The smallest interval system that contains all the possible crisp models is described by:
with
Suppose also that the true dynamics are contained within the interval model, with the following B vector
Then, based on simulated values of ω(k) and υ(k) and the input values of n d shown in Figure 5a, the application of the KF and IKF algorithms to the nominal model (Equations 5 to 7) and interval model (Equations 10 to 12), respectively, yield the estimates of the yaw angle shown in Figure 5b, in which the simulated true heading of the vessel is plotted as well.
It is clearly seen in Figure 5b how the nominal-system KF estimate lies within the IKF boundaries, and this holds for KF estimates obtained using any crisp model contained in the interval model. However, another phenomenon is evidenced: that of the rapid separation of the IKF boundaries. A consequence of the inclusion property of IA, which on the one hand is the raison d'être of most interval analysis applications and on the other constitutes one of the major practical difficulties in implementing IA-based algorithms, rendering results overly pessimistic and of little practical value. Hence, the consequently excessively conservative IKF boundaries may be attributed not to the theoretical framework of the algorithm laid out by Chen and co-researchers (Reference Chen, Wang and Shieh1997) (Reference Chen, Xie and Shieh1998), but with the way in which the IA calculations are implemented.
With regards to IA implementation, the IKF simulation was carried out with the aid of the open-source extension of MATLAB for IA, INTerval LABoratory (INTLAB), developed by Rump (Reference Rump1999). Numerous programming languages now contain libraries that extend the basic variable types to include intervals, and incorporate routines for carrying out IA operations (Kearfott, Reference Kearfott1996). Albeit in themselves highly efficient, the sharpness of the results obtained in a computation involving IA is highly dependent on the particular numerical formulation adopted, whereby even with the aid of such software, a naïve use of IA may provide results that are of no practical use. In fact, the cardinal reason for this reliance of the resulting interval on the precise numerical formulation adopted to compute some function is adequately known as the dependency problem of IA.
This phenomenon can be illustrated with a simple example: consider the interval x I=[−1,1]; the naïve computation of x I−x I yields [−1,1]−[−1,1]=[−2,2] whereas clearly x I−x I={y−y: y∈x I}=0. This failure of the cancellation law occurs because IA does not take into account that both intervals represent the range of the same variable, leading to highly overestimated bounds for the solution (in this case, infinitely overestimated). In higher dimensions, ie with vector operations, the dependency problem is exacerbated as what is known as the wrapping effect, which occurs when intermediate results consisting of sets of arbitrary shapes must be enclosed into boxes in order to be described by interval variables. In fact, Neumaier (Reference Neumaier1993) attributes this phenomenon to be the main difficulty of applying IA to dynamical system simulation.
Overestimation caused by the dependency effect occurs because each occurrence of a variable in a mathematical expression is implicitly assumed to be independent; therefore, this effect is suppressed if the expression can be reformulated so that every interval variable appears only once. Though this may not be possible for complicated expressions, one can in general minimise the occurrence of a single variable at a time in the expression. Consider, for example, the computation of the first element of the a priori error covariance matrix PI−(k)by the IKF algorithm for a second-order system:
where the shorthand notation pk≡p(k) has been used for convenience. Then, noting that x I(y I+z I)⊆x Iy I+y Iz I, for x Iy I, z I∈I(R), the expression in (14) can be reformulated in the following ways, each one taking advantage of a different factorization and yielding a different interval enclosure for p11,kI−:
the last three additionally making use of the symmetry of the error covariance matrix. It is not clear which one of these formulations, if any, yields the tightest enclosure, and even so, it would depend on the particular interval values taken on by each variable, which vary for each iteration. However, taking the intersection of all the intervals f i yields the narrowest enclosure of ${\rm p}_{11,{\rm k}}^{{\rm I} -} :\,\cap _{i = 1}^4\, f_i \subseteq f_i \forall i$. Using this approach for every component of each of the IKF equations for the Springer interval model (Equations 10 to 12), the improved IKF boundaries obtained are shown in Figure 6b. Also shown are the boundaries obtained using naïve IA, the KF estimate based on the nominal model (Equations 5 to 7), the simulated real trajectory, and the average of the improved IKF boundaries.
It is clear again from Figure 6b that without any special treatment, the widths of the interval bounds grow very quickly; however, in comparison to the improved boundaries it is evidenced that they encompass not just the actual set of possible KF estimates, but a much larger set that dwarfs the former, rendering the bounds meaningless in practice. The simple treatment carried out here to reduce the dependency effect has shown a significant decrease in the overestimation of the interval widths. Still further reduction may be obtained by applying more sophisticated techniques. For example, Alefeld and Mayer (Reference Alefeld and Mayer2000) show that a centred first-order Taylor representation of a function provides a quadratic order of approximation of the actual range, while Neumaier (Reference Neumaier2009) discusses exploiting monotonicity properties. Another strategy may consist of splitting the input intervals into smaller subintervals, obtaining an interval-enclosure of each one, and taking the union of these, since it has been shown to reduce the dependency effect (Daumas et al., Reference Daumas, Lester and Muñoz2009). Techniques have also been proposed to reduce the aforementioned wrapping effect of IA, such as Neumaier's (Reference Neumaier1993) use of ellipsoids to compute affine transformations of vectors in which both the vectors and the transformation matrix contain interval elements, or Kühn's (Reference Kühn1998) use of zonotopes, a kind of polytope, to approximate the system states.
Of particular importance is the computation of a good enclosure for the inverse interval matrix for the Kalman gain: since an interval matrix containing any singular point-valued matrix is not invertible, overestimation of intervals may result in the computed matrix being non-invertible, preventing further application of the IKF algorithm. Owing to this, several workaround strategies to bypass this difficulty have been proposed. For instance, Chen et al. (Reference Chen, Wang and Shieh1997) proposed using a “suboptimal IKF” in which the inverse interval matrix is replaced by its worst-case inversion, which is an ordinary (non-interval) matrix. More recently, Xiong et al. (Reference Xiong, Jauberthie and Travé-Massuyès2012) have shown how to use the set inverter via interval analysis (SIVIA) method, an algorithm used to search the pre-image of a given set under a given function (Jaulin and Walter, Reference Jaulin and Walter1993), to provide a tighter enclosure for the interval Kalman gain, while maintaining guaranteed bounds. In the wider field of interval analysis, efficient computational schemes for interval matrix inversion are being investigated, with work having recently been carried by Nirmala et al. (Reference Nirmala, Datta, Kushwaha and Ganesan2011).
Recently, in order to address the conservatism inherent to IKF, Ahn et al. (Reference Ahn, Kim and Chen2012) have derived a new filtering scheme from the IKF in which the interval uncertainties in the system model are incorporated into the covariances of the statistical model. Though this new filtering scheme is able to take into account interval model uncertainty, it however provides a point-valued estimate rather than the rigorous bounds of the IKF.
IA is also being applied to other robust filtering techniques. Kieffer et al. (Reference Kieffer, Jaulin and Walter1998) presented a state estimation technique in which system and measurement uncertainty are described as random sequences but bounded by intervals, rather than by the traditional stochastic modelling assumptions of KFs. This new technique, based on IA and set inversion, returns as estimate at each time-step a set guaranteed to enclose all values of the state that are consistent with the observations. Prediction and correction phases are alternated in a way reminiscent of the KF. Whereas Reece (Reference Reece2000) proposed a new filter, the Biscay distribution filter, which combines IA with statistical KF estimation methods, Ashokaraj et al. (Reference Ashokaraj, Tsourdos, Silson and White2004) demonstrated the combined use of IA and the EKF for robot navigation, in which localisation estimates from the two methods were fused to create a more robust estimator.
Soon after Chen et al. (Reference Chen, Wang and Shieh1997) published their IKF algorithm, an extended version for non-linear systems, the extended IKF (EIKF), was derived, and its usage was demonstrated in a missile-tracking simulation problem (Siouris et al., Reference Siouris, Chen and Wang1997). He and Vik (Reference He and Vik1999) extended the study of the EIKF by applying it to a simulation of an integrated GPS-INS system for aircraft navigation. Additionally, the simulation of a ship navigation system using the IKF to fuse GPS readings with dead reckoning calculated from compass and velocity measurements was reported by Tiano et al. (Reference Tiano, Zirilli, Cuneo and Pagnan2005).
Though the IKF computes the bounds of the KF estimates of every point-valued system model contained in the interval model, in practice a unique value of the estimate is often needed as well. Chen et al. (Reference Chen, Wang and Shieh1997) had suggested using a weighted average of the boundaries, or more simply, the average of the boundaries. The average of the IKF boundaries for the Springer yaw estimation is shown in Figure 6b, where it can be seen to roughly coincide with the KF estimate obtained using the nominal system model. In another study, a fuzzy logic inference scheme was applied to the IKF (Chen et al. Reference Chen, Xie and Shieh1998) to deliver point-valued estimates. Later on, Weng et al. (Reference Weng, Chen, Shieh and Larsson2000) went on to describe the use of evolutionary programming techniques on an IKF to firstly, reduce the interval estimate at each iteration to an “optimal interval”, and secondly, find a nominal value within that interval which best represents the state estimate for practical purposes.
6. CONCLUDING REMARKS
While most USVs have typically been remotely operated, the demand for USVs with increasing levels of autonomy is on the increase. With the large number of emerging USV platforms and the wide use of the KF in navigation and for sensor fusion in general, research is currently being undertaken in developing robust navigation techniques for USVs, as changing maritime conditions, payload, etc, characterize these vessels. These translate into uncertainties in the dynamical model that can usually be easily bounded by intervals. The IKF is a natural extension of the KF for interval system models that provides rigorous bounds to the estimate in the face of such bounded model uncertainty.
In this paper, the IKF and KF algorithms are applied to the yaw angle estimation of the Springer USV, revealing the inclusive nature of the IKF estimate and its overly pessimistic yield when naïve IA is used. Undesirable interval widening is a common problem in all applications that require IA. Techniques are being developed to reduce this effect in the broader field of interval analysis to solve problems that require rigorous computing, where these may be adapted for use on the IKF algorithm. The simulation also shows how the IKF may substitute the regular KF algorithm as point estimates can be inferred from the boundaries if needed. The paper also registers the limited use the IKF has had since its inception, even though several improvements have been proposed and studied through computer simulations and its usage as a state estimator for systems with bounded model uncertainty being regarded favourably despite the difficulties of its implementation. It is intended that this paper will provide insight into the IKF's purpose and applicability, as well as presenting the reader with a broad picture of the current state of the art regarding improvements to the IKF developed, with proposals for further investigation. Some of these are focused on improving boundary sharpness and others on inferring reliable point-valued estimates from the intervals. The inclusion of IA libraries in programming languages makes it possible to implement these strategies. The versatile payload of many USVs allows for full-scale computers to be used on board, as well as being able to host a wide variety of sensors, making USVs ideal test-beds on which to assess the viability of implementing IKF-based solutions in real systems, providing the incentive for its widespread use, while at the same time addressing the necessity of the USV market for increasingly robust navigation systems.
APPENDIX
Consider the following linear stochastic-deterministic system
where x(k), u(k), y(k) are vectors of the adequate dimensions that represent, respectively, the system state, the (deterministic) system input, and the measurement vector at time-step k, and ω(k) and υ(k) are white noise sequences with zero-mean Gaussian distributions with known covariances cov(ω)=Q, cov(υ)= R, and that E[ω(l) υT(k)]=0l,k, E[x(0) ωT(k)]=0, E[x(0) υT(k)]=0∀k.
Then the KF equations provide a statistically optimal (unbiased and minimum error variance) estimate of the true state vector, and can be written as a set of recursive equations (Figure A1), from which the state estimate at each time-step is obtained from the previous estimate and the new observed measurement, assuming initial estimates for x(0) and PI(0).