APPLICATION OF THE KALMAN FILTER TO ICEBERG MOTION FORECASTING Christophe Simon B. Sc., Ecole Superieure d'Ingenieurs de Marseille, 1988 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 1990 © Christophe Simon In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada Department Date DE-6 (2/88) Abstract The objective of this study is to develop an application of the Kalman filter for filtering and forecasting iceberg positions and velocities in order to calculate the risk of impact against a fixed structure or stationary vessel. Existing physical and probabilistic models are reviewed. Physical models are essen tially based on the response of the iceberg to the forces acting on it. Statistical models forecast the motion of the iceberg based on past observations of the trajectory. A probabilistic iceberg trajectory model is used in this study so that uncer tainties in the trajectory forecast can be explicitly included. The technique of Kalman filtering is described and applied to forecast future positions and velocities of an ice fea ture, including uncertainties.The trajectory forecast combined with a risk calculation, yields the probability of impact and the probability distribution of the time until impact. Numerical simulations are provided to demonstrate and validate the procedure. ii Table of Contents Abstract ii List of figures v Nomenclature vi Acknowledgements vii1 Introduction 1 2 Existing iceberg trajectory forecasting models 4 2.1 Physical models 4 2.1.1 Dynamic model 5 2.1.2 Kinematic model 6 2.2 Shortcomings of present approaches 7 2.2.1 Physical input uncertainties2.2.2 Collection of field data 8 2.3 A probabilistic approach 9 3 Iceberg motion filtering and forecasting 13 3.1 A brief introduction to the filtering problem3.2 Interest for iceberg motion forecasting 14 3.3 A recursive approach of data filtering: the Kalman filter 15 4 Application of the Kalman filter to iceberg motion 19 iii 4.1 Formulation of the Kalman filter 19 4.2 Choice of the process vector dimension 23 4.3 Trajectory filtering 25 4.4 Noise calculation 6 4.5 Trajectory forecasting 9 5 Risk of impact calculation 32 5.1 A numerical algorithm5.2 Application of the filter to a real trajectory 35 5.3 An example of risk calculation 42 5.4 Comparison of risk calculation and simulation 50 6 Conclusion 56 Bibliography 8 Appendix A : Flowchart of the filtering software 61 Appendix B : Flowchart of the forecasting software 62 Appendix C : Flowchart of the risk calculation software 63 iv List of Figures 3.1 Algorithm of the Kalman filter 18 4.2 Computed trajectories 24 4.3 Convergence of errors during filtering 27 4.4 Real iceberg trajectories 30 5.5 Actual iceberg trajectory used in example 36 5.6 One-step-ahead forecasted position 7 5.7 The newly observed position 38 5.8 Filtered position 39 5.9 Filtered trajectory 40 5.10 Forecasted trajectory 1 5.11 Increase of uncertainty perimeters (r.m.s errors) . 43 5.12 Increase of position uncertainty with time 44 5.13 A computed trajectory 45 5.14 Corresponding "observed" trajectory 47 5.15 Trajectory forecasting and position error perimeters 48 5.16 Forecasted positions and corresponding times 49 5.17 Calculated risk distribution 51 5.18 Risk distribution obtained with a simple simulation! 53 5.19 Uncertainty over the direction of the future trajectory 54 5.20 Risk distribution obtained with a complete simulation! 55 v Nomenclature [18] = reference to item 18 in Bibliography a = constant in the autocorrelation function Ak = matrix relating two steps of the state vector in absence of the forcing function at time tk d — distance between center of iceberg and center of platform Ck = observation matrix connecting measurement and state vectors at time tk E[] = expected value operator ER = mean-square error ejt = error vector between forecasted and actual observations at time tk FD = fluid drag force K~k = gain of the Kalman filter at time tk ko = frictional drag coefficient ki = real coefficients in kinematic model n(T) = number of impacts between time=0 and time=T p[z = Z) = probability that the random variable z has value of Z p(z = Z\y = Y) = probability that the random variable z has value of Z given that the random variable y has value Y Pk = error covariance matrix at time tk P£ = error covariance matrix at time tk for one step forecasting Qk = covariance matrix of to* i?() = autocorrelation function of iceberg velocities vi Rk = covariance matrix of Vk i(t) = vector position of the iceberg at time t Rot() — a rotation matrix function Sj = influence area of an iceberg Sc = cross-sectional area of an iceberg s,t = time constants tk = discrete time u(t) = velocity of the iceberg at time t u(t) = best estimate of the velocity at time t U{ = discrete values of the velocity Vk = observation white noise vector Wk = process white noise vector x(t) — signal at time t x(t) = filtered value of the signal at time t Xk = process vector at time tk xk = best estimate of the process vector at time tk xl~ = best estimate of the process vector knowing the process prior to time tk Z(k) = observation vector Z(k) = estimate of the observation vector a, 9 = angles 7 = time constant of the autocorrelation function St = small time increment if - angle p = density of fluid (i — coefficients vii Acknowledgements I wish to thank Dr. A. B. Dunwoody, my supervisor, for his guidance and support throughout the course of this research. I also must thank Mona El-Tahan of C-CORE who supplied the iceberg trajectory data, and the Natural Sciences and Engineering Research Council who sponsored the work through an operating grant. Finally, I wish to thank Susanne for her help in editing this thesis. viii Chapter 1 Introduction Icebergs are a major threat for offshore structures and boats in eastern Canada, and especially in the Northwest Atlantic Ocean and in the Labrador Sea. They have been a subject of interest as a danger for navigation since the last century. One can think about the Titanic disaster in 1912. To prevent offshore accidents the International Ice Patrol collects information con cerning drifting ice from boats, aircraft, satellites and other sources. Iceberg trajectories are then set up to provide the consequent warnings. Even though iceberg mapping provides a good measure of the risk encountered in terms of ice density and short term motion forecasting, it appears not to be sufficient if applied to fixed offshore structures. The renewed interest in the detailed movements and forecasting of icebergs began with the oil and gas exploration in the continental shelf off eastern Canada. The exploration activity on the Grand Banks and in the Labrador Sea has required more detailed tracking and forecasting of iceberg movement in the vicinity of drilling activities. The predominant type of vessel for oil exploration in the area is the semi-submersible drill rig. This type of platform cannot withstand an impact from an iceberg. Even fixed oil-production structures are not able to withstand impact from the largest icebergs. To minimize the threat posed to stationary floating and fixed facilities, it is necessary to forecast the trajectories of all icebergs in the vicinity of a facility over periods from 10 to 24 hours into the future. Appropriate actions have to be taken if it appears that an 1 Chapter 1. Introduction 2 iceberg might hit. The emergency procedures vary depending on the type of operation and mooring. It can take several hours before the platform can be released and the well secured. In some circumstances, it might be preferable to tow the iceberg away. In the worst case, it would be necessary to evacuate personnel by boats and helicopters. It is obvious that the operators of the platform do not wish to initiate any of these safeguards if they are not necessary. Therefore, it is important to be able to evaluate the risk that an offshore facility will be hit by an iceberg which is in the vicinity before the appropriate safeguards are implemented. It should also be easy to update the level of risk with any new development of the iceberg's motion. Currently, any oil exploration activity must include an iceberg monitoring program. Iceberg motions are tracked by radar from the rig. The experience and judgement of the iceberg monitoring personnel are then used to determine the risk of an impact. In order to help the iceberg monitoring personnel, numerous efforts have been made in the last 20 years to predict the trajectories of icebergs based on the physical or statistical modelling of iceberg trajectories. Two of these models also explicitly evaluate the risk of an impact by incorporating the growing uncertainty in position and velocity of the iceberg with forecasting time. The existing trajectory models will be reviewed before focussing on one model for detailed development. The trajectory model to be described uses past positions of an iceberg to predict the future trajectory, and hence its risk of impact against the facility being protected. Kalman filtering is an established technique for predicting future behaviour of a stochastic system from a history of observation, including noise in the measurement of the process as well as the "process noise" driving the system. The Kalman filter will be applied to trajectory forecasting of icebergs. Chapter 1. Introduction 3 Of course, the ultimate goal of any iceberg monitoring program must be the evaluation of the risk of impact. The risk equation derived by B. Dunwoody [11] uses the forecasted trajectory and the uncertainties on positions and velocities given by the Kalman filter to evaluate the risk of impact. Finally the results of a numerical simulation will be described which illustrates and validates the method. Chapter 2 Existing iceberg trajectory forecasting models Since the early seventies, several types of prediction procedures have been developed. They can be classified into two main categories: 1. evaluating the influence of the physical environment on the iceberg, and deducing its trajectory step by step as a function of current and wind fields, 2. auto-regressive models using the known trajectory as the main source of input. The first type of method has been the only one used until the beginning of the eighties. Auto-regressive models are more recent. They have been introduced to handle the large uncertainties existing in the wind and flow fields, as well as in the physical charac teristics of icebergs. 2.1 Physical models A physical model uses winds and currents to calculate iceberg motions either with a dynamic or a kinematic approach. It requires certain informations about the iceberg and its surroundings. Environmental data can be obtained from different sources. The wind is measured on the platform or ship. Forecasted data are provided by the Canadian Atmospheric Environment Service or other sources. Current meters are deployed from the platform. They provide information concerning the flow field in the vicinity of the structure, sometimes at different depths. 4 Chapter 2. Existing iceberg trajectory forecasting models 5 Shape and mass of icebergs can be reasonably estimated from visual observation. Extensive data have been collected for this purpose in the Labrador Sea [18]. They relate aircraft observations to total mass, dry and wet areas. Physical parameters such as bottom topography and water depth are collected during site oceanographic surveys. 2.1.1 Dynamic model The earlier studies have mainly used the dynamic equation to describe iceberg trajecto ries. Since then, several procedures have been developed [10,14,18,23]. The iceberg is subjected to several forces. Those forces vary in amplitude and direction depending on the current and wind. The dynamic equation expresses that: Mass of Iceberg x acceleration = Sum of external forces. A model based on this equation has to take into account the horizontal components of the forces due to water drag, wind drag, Coriolis acceleration and sea surface slope or pressure gradient. The fluid drag force on the iceberg can have two different forms : where p is the density of fluid, kr> is the frictional drag coefficient, u is the velocity and Sj is the influence area. This equation applies when the frictional drag due to the viscous stress acting at the surface of the iceberg, is predominant. However, when the form drag is predominant, the equation that applies best is FD - pkDSju (2.1) FD =^pCDScu2 (2.2) Chapter 2. Existing iceberg trajectory forecasting models 6 where CD is the form drag coefficient and Sc is the cross-sectional area. The main components of the water movement are the geostrophic current, the wind-driven current, the inertia current and the tidal current. Their magnitudes and directions vary with time and location. The water column usually has to be divided into several layers to take into account the variability of water current with depth [18,23]. The water drag on an iceberg seems to be often the most important forcing term, but the wind drag has sometimes to be considered, depending on the shape of the iceberg and on atmospheric conditions. A wind greater than 30 knots and blowing from a constant direction can be a major factor. However, in general, the ratio of iceberg velocity relative to wind speed is only 0.03 [20]. The Coriolis force acts on the iceberg and on the water in which the iceberg floats. Therefore, the Coriolis force created by the movement of the iceberg is balanced by the pressure gradient in a geostrophic current. The pressure-gradient force for each layer can be calculated as the Coriolis force of the current, while the velocity of the iceberg relative to the water gives the other part of the overall Coriolis force. One has also to consider the inertia drag due to the acceleration of the object in terms of a hydrodynamic mass to be added to the actual mass of the iceberg. 2.1.2 Kinematic model As opposed to the dynamic point of view, the kinematic model expresses that the iceberg velocity is a function of the velocities of the different influences [7]. This type of model has some success for two reasons: 1. it appears from observations in the east cost of Canada that icebergs generally respond quickly to external forces (time scale of about one hour or less), provided that the considered influence is significant [7], Chapter 2. Existing iceberg trajectory forecasting models 7 2. and there is a good correlation between the motion of water and the drift of an iceberg, which has 7/8 of its volume below water. The vector velocity of the iceberg u(t) is empirically calculated as a function of wind and current velocities: u(t) = ki Rot((pi) (2-3) where u^) represents the vector velocity of the ith influence, Rot is a rotation matrix, is the phase angle and fc, is the empirically determined coefficient [7]. The water influence is divided into wind generated current, tidal current and so forth. The values of the empirical coefficient k and the phase angle (p are related both to the size of the iceberg and to the rate of change in the external influence. The smaller the iceberg or the rate of change, the closer is k to 1 and if to zero. 2.2 Shortcomings of present approaches None of the existing physical models has been extensively used in the field, because of practical difficulties. 2.2.1 Physical input uncertainties The accuracy of the dynamic model is dependant on the accuracy of the physical in put. Equations of motion remain sensitive to errors in drag coefficient and other shape parameters, as well as towards errors in the current and wind values. The external shape of an iceberg can tell about its characteristics, but would only give a rough approximation of the reality. Chapter 2. Existing iceberg trajectory forecasting models 8 Current and wind values measured in the vicinity of the platform are biased by its presence. The same type of accuracy problem applies with the prediction of wind and current. Even if long term current, mean flow and tides are predictable, one cannot forecast with sufficient accuracy local perturbations and low frequency water motion. Atmospheric parameters may also be highly spacially and temporally uncorrelated. Furthermore, the low pressure generated currents can be more significant in affecting iceberg drift than wind or wind-generated currents [7]. Different approaches have been proposed to compensate for inexact values. Dempster [8] suggested adjusting the drag coefficient of the iceberg in order to obtain the best fit between the calculated and the observed trajectories. He used a mix of dynamic and kinematic considerations. Smith and Banke [22] used an hincast model of motion, and adjusted the air and water drag coefficients to compensate for imprecise knowledge in the physical characteristics of the iceberg. 2.2.2 Collection of field data For practical applications, dynamic and kinematic models are of no use unless the input parameters are available. The collection of field data presents several complications: 1. equipment mounted on the platform are subjected to failure or environmental re strictions such as storms or pack ice, 2. information collected on the platform is valid only in the vicinity of the structure. Remotely based strings of current meters are a possible but very expensive solution for this problem [2]. Chapter 2. Existing iceberg trajectory forecasting models 9 If wind and current are not avalable, the only possible source of information is the past trajectory of the iceberg. One can try to deduce wind and current from known radar positions, and then apply the obtained values to forecast the future motion [14,15]. Here again, existing models suffer from uncertainties and variability in data input. Dynamic and kinematic procedures tend to use exact and difficult procedures with inexact data. This level of accuracy seems inappropriate to the forecasting problem, because of uncertainties in wind and flow fields, poor knowledge of icebergs shapes, unreliability of current and wind forecasting. Exact models are certainly more related to scientific interest than to field application. Furthermore, the physical models previously described don't allow for the uncertain ties for position and speed necessary for risk calculation. There lies an alternative to dynamic and kinematic methods in the following approach. 2.3 A probabilistic approach As mentioned previously, another and more recent field of interest is to compute iceberg trajectories in a probabilistic manner. The idea is that if we can't really know where the iceberg will go, we should be able to measure our error, or the area where the iceberg has the most chances to go. Using the same approach, C.Garret [12] writes the "best estimate" of the future velocity u(t) in terms of the known past values the coefficient £„ being deduced from a set of equations expressing the minimization of the mean-square error: V,(t) = E„Cn«, n (2.4) ER = (u{t) - u{t)Y (2.5) Chapter 2. Existing iceberg trajectory forecasting models 10 giving the set of equations: £m "n um (m = u(t)un (2.6) The solution can be found given the covariance matrix un um and the covariance vector u(t) un. In order to solve the equation of the mean-square error, one has to know in what ways the velocities are related to each other. The autocorrelation of u is defined as: R{T) = u(t)u(t + T)/u2 (2.7) Therefore, rt+tn (2.6) « Sm R(tm - tn) (M = / R(r)dr Garrett takes an exponential autocorrelation function expressing a stationary Gauss-Markov process: R(T) = e -7T (2-8) The exponential form indicates that the velocities become less and less correlated as the time separation between samples increases. Because of noise in the velocity data, the autocorrelation behaves more like (2.9) R(r) = ae-<T if r > 0 R(0) = 1 The constant a expresses the presence of noise in the velocity data. Because a < 1, the value of the autocorrelation function is reduced. Chapter 2. Existing iceberg trajectory forecasting models 11 a and 7 are calculated using known trajectories and field data (wind and current) processed through lagrangian and eulerian analysis. Data had been collected during an oceanographic survey of the area. Solving equation (2.6) gives the set of coefficient £„. The mean square position error normalized with mean square velocity is: ER(t) = (x(t) - x(t)Y/u* (2.10) and is calculated using a and the coefficients („. Error calculation is an important step towards risk prediction. Garrett supposes that the probability of the iceberg to be at a distance d from the forecasted position is gaussian, ie can be calculated using the mean square position error. The probability of impact is then deduced using a simple procedure based on: • the probability of presence in a circle around the platform, • the expected transit time and the average passage time through the circle. The probability of speed is taken independent of initial velocity perturbation. B. Dunwoody [11] takes into his forecasting scheme a Wiener process to express the rate of increase of uncertainty in velocity with time. This simple model has a variance increasing linearly with time. The expected velocities and positions are deduced from the last known position. Risk calculation is based on the expression of the joint probability distribution of position and speed. Positions and speeds are then integrated around the wellsite, giving risk of impact and distribution of speed upon impact. Once again, the error relative to the expected position is at the base of risk calculation. After having reviewed the existing prediction models, one can already conclude that Chapter 2. Existing iceberg trajectory forecasting models 12 an exact calculation of the movement of an iceberg cannot be carried out in the field with accuracy. Wind and current are usually not available from remote sites in real time. Except for mean values, they are spacially and temporally uncorrelated, and therefore unpredictable. It is possible to get a probabilistic image of the perturbations in flow and wind fields for a given area. Other sources of information concern the close vicinity of the platform (local current and wind), as well as radar tracks of icebergs. Until recently, the main field of research has been the physical description of the ice berg's trajectories. A more realistic approach of the iceberg motion forecasting problem would be to develop a probabilistic model. This model could use the past trajectory as the only input, in order to give a forecasting as accurate as possible in some sense. It would also give the corresponding uncertainties over positions and velocities. Chapter 3 Iceberg motion filtering and forecasting The auto-regressive models do not give significantly poorer forecasts in practice than dynamic or kinematic models. They have two main advantages: 1. they require less data to implement, which is of a major importance for practical use, 2. they can incorporate uncertainty in the forecast. It is therefore reasonable to develop an auto-regressive model to generate trajectory forecasts, and to couple it with Dunwoody's risk calculation. This implies that the forecasting would also generate the first and second moments of position and velocity. Both forecasting and uncertainties calculations can be done by "filtering" a sequence of measured positions of an iceberg. 3.1 A brief introduction to the filtering problem In signal theory, the traditional purpose of a filter is to separate the signal from the noise. In communication, a filter gives a certain frequency response.' In the forties, Norbert Wiener considered a different type of problem, where signal and noise are both noiselike in character, with possible overlap in their spectra [25]. In this case, perfect separation is not possible. One can only expect a "best" estimated signal, according to the chosen performance criterion. Wiener's assumptions were the following: 13 Chapter 3. Iceberg motion filtering and forecasting 14 1. The filter is linear and stationary with time, 2. the filter input is an addition of signal and noise, whose covariances are stationary with time, and with known auto- and cross-correlation functions, 3. the output's covariance is stationary with time, 4. if t is the present time, x(t) is the signal and x(t) the filtered value, the performance criterion is the minimization of the mean square error: ER((x(t) - x(t + s))2) If s — 0, this is the filtering problem; s > 0 is the prediction problem, and s < 0 is the smoothing problem. Garret's prediction scheme is indeed a direct application of the so-called Wiener filter. This is the weighting function approach. The auto-correlation of the iceberg's velocity has been calculated using extensive field measurement . 3.2 Interest for iceberg motion forecasting The iceberg trajectory forecasting problem can be treated as a form of filtering problem. The "signal" to be filtered is the trajectory of the iceberg. The "noise" from which the signal must be extracted is 1. one cannot measure iceberg speed directly, 2. errors exist in position determination. The dynamic system of the iceberg motion can be expressed as a difference equation with "noise" input being fluctuations in the environmental loads imposed on the iceberg. Chapter 3. Iceberg motion altering and forecasting 15 Known wind and current mean values produce a deterministic "noiseless" motion. All other unknown or non-measurable influences such as turbulence, high frequency wind and water motion cause "noisy" behaviour of the iceberg. Of course, the Hmit between mean and noisy influences is not obvious, but deterministic and noisy parts can either be treated separately or mixed together to produce a different auto-correlation function. The Kalman filter is an efficient tool for the filtering of a noisy signal. 3.3 A recursive approach of data filtering: the Kalman filter Following Wiener's work, a major contribution in the filtering problem was made by R.E Kalman in 1960 [16]. The Wiener filter gives the "best" of present and future process values using the weighting of all known past measurements, assuming that the auto-correlation of the process is known. Kalman kept the same performance criterion, ie the minimum mean square error, but introduced a recursive technique for data processing. In other words, the Kalman filter is self-adaptive to new input. All past values are step-by-step included in the filter's arithmetic. This is the main difference between the Kalman filter and the weighting function approach. The recursive technique needs both process and measurement procedures to be mod eled in vector forms. The following description of the Kalman filter refers to a discrete time formulation, which is our concern as we will see later. We assume that the random process can be modelled as: and the measurement process as: = Akxk + wk (3.11) Chapter 3. Iceberg motion filtering and forecasting 16 zk = Ckxk + vk (3.12) where • xk = (nxl) process state vector at time tk • Ak = (nxn) matrix relating two steps of the state vector in absence of the forcing function • wk = (nxl) white noise vector with known covariance structure • zk = (pxl) observation vector at time tk • Ck = (pxn) observation matrix connecting measurement and state vector • vk = (p x 1) measurement vector error, white noise with known covariance structure, and uncorrelated with wk • Qk and Rk are the covariance matrices of wk and vk: Qk = E(wkwl) Rk = E(vkvl) where T denotes the transpose. We assume that: 1. E(wkwJ) = 0 for k ^ / 2. E{vkvJ) = 0 for k ^ I 3. E{wkvJ) = 0 for all k and /. Chapter 3. Iceberg motion filtering and forecasting 17 Let Xk be the "best" estimate of the state vector at time tk, and x]^ the "best" estimate of the state vector knowing the process prior to time tk-The idea is to express Xk with the following form: Xk = xl + Kk(zk - CkXk ) (3.13) Let ejt = zk - CkXk (3-14) be the error between the forecast of the present output CkX^ and the present measure ment Zk- ek is weighted with the gain matrix K~k-In the Kalman filter theory, the gain of the filter is calculated in order to minimize the variances of the estimation error for the state vector, the terms in the major diagonal of: Pk = E[(xk - Xk)(xk - xk)T] (3.15) Another important matrix is: Pk~ =E[(xk-Xk)(xk-Xk)T] (3.16) whose elements are the variances of the forecasting error at time r^. A complete descrip tion of the Kalman filter theory can be found in [1,5,16,17]. The detailed algorithm of the filter is given in figure (3.1). Chapter 3. Iceberg motion altering and forecasting Previous estimate Error covariance A Prediction c Estimate of present output ik = Ckxk Predicted covariance Error Present output Error covariance for updated estimate Pk = (I- KkCk)Pk-Current, estimate xk = x\~ + Kkek Correction A'jtejfe Kk = Filter gain ircI(ckFTCl + Rky Figure 3.1: Algorithm of the Kalman filter Chapter 4 Application of the Kalman filter to iceberg motion 4.1 Formulation of the Kalman filter For iceberg motions, we are concerned with a discrete time model, since position mea surements (radar plot) are taken on a discrete time basis. The state vector xk includes position, speed and possibly acceleration of the iceberg. Therefore, it can be either a (4 x 1) or a (6 x 1) vector. The horizontal plane is (x, y): xk - , . dx dy where x = — and y = — dt y dt xk Vk \ y*) For a discrete time model, the new position is the old position plus old velocity times time increment. The new velocity is the old velocity plus the integral of acceleration. The process equation is therefore Xk+i = Xk + Txk ik+i = ik + w-i J/fc+i = Vk + Tyk yk+i = yk + W2 (4.17) 19 Chapter 4. Application of the Kalman filter to iceberg motion 20 where T is the time interval between two observations, or in vectorial form: (the subscript h of Ak is omitted if A is constant) (4.18) with A = and Wk ( 1 T 0 0 ^ 0 10 0 0 0 1 T 0 0 0 1 The process noise covariance matrix is therefore V ' 0 ^ 0 Qk = E[wwT) = 0 0 0 0 0 E[w\] 0 0 0 0 0 0 ^ 0 0 0 E[wl) ) The observation equation is and V y J Vk V J + \ v* } (4.19) (4.20) Rk = E[vvT] = E[v\) 0 0 E[vj] (4.21) Chapter 4. AppHcation of the Kalman filter to iceberg motion 21 If acceleration is included, the state vector becomes: I r ^ Xk Xk Xk Vk yk \ y« i i dx . ,. dy where x = — and y = — dt * dt The process equation is therefore Xk+i — Xk + Tik + ^-xk ifc+i = xk + Txk Xk+l = Xk + W! Vk+i = yk + Tyk + %yk Vk+i = yk + Tyk Vk+i = yk + w2 (4.22) or in vectorial form: *fc+i = Axk + Wk (4.23) Chapter 4. Application of the Kalman filter to iceberg motion 22 with A — ( 1 T f 0 0 0 \ 0 1 T 0 0 0 0 0 1 0 0 0 0 0 0 1 T ^ 0 0 0 0 1 T y 0 0 0 0 0 1 The process noise covariance matrix is now: and Wk = Qk = E[wwT] = I 0 ^ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E[W\) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E[w (4.24) The observation equation is y J 1 0 0 0 0 oN 0 0 0 1 0 0 , / \ xk Xk Xk Vk Vk + \v2 J (4.25) and v and R remain unchanged. Chapter 4. AppHcation of the Kalman filter to iceberg motion 23 The proposed discrete time physical model is the simplest available. More sophisti cated models (eg including known currents, etc) could be substituted if warranted. 4.2 Choice of the process vector dimension In the process equation, the noise can be applied either on the velocity or the acceleration of the iceberg. The better model is the one giving results as close as possible to reality. Random iceberg trajectories have been created to compare with real trajectories. The artificial trajectories are calculated step-by-step from a given initial position Xo, using the process equation: xk+T. = Axk + w The noise w is taken randomly in a normal distribution with zero mean and given (22[iWi], -Efu^]) variance. The time interval between two steps is one hour, as this is the most common value for actual radar tracking data. Two sets of random trajectories have been created. In the first set, the noise has been applied on velocity. In the second set, trajectories have been calculated with noisy acceleration. In both sets, noise variances have been chosen in such a way that the distance between two positions remains realistic. A typical output is presented for each set in figure (4.2). Comparison with real trajectories shows that a noisy velocity is the more realistic of the two. For a time scale of one hour, the rate of change of position is too slow with a noise applied on acceleration. In this case, a realistic behaviour would mean a very large value for the noise. The velocity model is therefore kept to proceed to the numerical integration of the risk in chapter 5. Chapter 4. AppHcation of the Kalman filter to iceberg motion 24 Figure 4.2: Computed trajectories Chapter 4. Application of the Kalman Alter to iceberg motion 25 4.3 Trajectory filtering Filtering iceberg trajectories means extracting the "best" estimate of the real trajectory from the measurement values. To implement the filter, estimation of the state vector JC,nit and covariance matrix Pinit are needed. In case of dimension 2, let (x(l), y(l)) and (x(2), y(2)) be the measured positions at time t\ and t2. An estimation of the state vector at time t2 can be deducted: y(2) \ l(y(2)-y(l)) ; The initial error covariance matrix is x(2) i(x(2)-x(l)) (4.26) P2 = E[(x2 - x2){x2 - x2)T] (4.27) Because x2 — x2 = -«i(2) Wl(l) - ±(Vl(2) - Vl(l)) -v2(2) \ w2(l) - ±(v2(2) - v2(l)) j (4.28) Let = < E\w\) = a E[w\] = b E[v\] = c E\vl) = d Chapter 4. Application of the Kalman filter to iceberg motion 26 a 0 0 \ a T 0 then P2 = (4.29) 0 0 b T \ 0 0 b 26 , J / x2 and P2 are the "best" estimates we can provide to the filter at time t2. At each successive time step, the filter uses only the measured position to calculate process values and errors, as shown in figure (3.1). As more and more measurement values are given to the Kalman filter, its errors and gain converge towards stationary values. Figure (4.3) compares the convergence of errors and gain related to x, either with the matrix P2 calculated above, or with P2 — (0). We obtain the same values after a few steps. In dimension 3, three measurements are required in order to calculate a first estimate of the process vector, including position, speed and acceleration of the iceberg. The convergence of errors and gain is obtained after 3 to 5 steps. 4.4 Noise calculation The Kalman algorithm requires the values of both process and measurement noise. Measurement noise depends on the accuracy of the radar. Typical values are 1 % of the range, and ±1 deg over the bearing angle [4]. The evaluation of the process noise is a more complicated problem. It covers all non-deterministic parameters influencing the motion of the iceberg. One method of obtaining this information would be to measure a set of iceberg trajectories,as well as wind and current values for the certain area and at a given period of time. This data can then be used to deduce the corresponding value of the noise. This is a laborious operation which Chapter 4. Application of the Kalman filter to iceberg motion 27 0.045 -r 2 4 6 8 10 12 14 16 time (hours) • initial conditions=(0) O calculated initial conditions Figure 4.3: Convergence of errors during filtering Chapter 4. Application of the Kalman filter to iceberg motion 28 requires a detailed survey of the area. Different icebergs can have different behaviours when driven by the same winds and currents. Therefore, it is desirable to use the motion of the observed icebergs as much as possible as the only source of information. In dimension 2, the process noise is appHed on velocity: ik+i = ik + wi Vk+i = yk + u>2 (4.30) Wi = T x ( acceleration on x) Therefore, 1 iy2 = T x ( acceleration on y) and E[ww2] is the variance of T times acceleration. To calculate the noise, the known trajectory can be used without treatment, by calculating velocity and acceleration versus time. The filter can also be used a first time with an arbitrary value for the noise. Then the noise over the filtered is calculated. The procedure is repeated until the value of the noise converges towards a final value. In dimension 2, a first estimation of the process noise using the observed trajectory gives better results. This appHes in the common case when the measurement noise is smaller than the process noise. Figure (4.4) shows three different iceberg trajectories observed off the east coast of Canada [24]. Locations were recorded at hourly intervals. The foUowing table gives for each of the trajectories: Chapter 4. Application of the Kalman filter to iceberg motion 29 • the noise calculated with respect to x i£[u>i], • the noise with respect to y JSfifj], • the overall noise \jE\w\)2 + E[u>l]2. Units are km2 h~4. noise on x noise on y overall noise Traj.l 0.07 0.07 0.10 Traj.2 0.20 0.08 0.21 Traj.3 0.06 0.11 0.13 The process noise is a good image of the "noisy" behaviour of the iceberg, and there fore of the predictability of its motion. The first trajectory presents sudden changes in both x and y directions, but on a small scale. The motion of the second iceberg is obvi ously more unpredictable on the x direction. The third trajectory presents more changes in acceleration in the y direction. The same methods of process noise calculation can be applied in dimension 3, where < Vk+i = Vk + w2 4.5 Trajectory forecasting The most important feature of the Kalman filter in developing the risk of impact calcu lation, is trajectory forecasting. At the end of the known trajectory, we obtain the "best" estimate of the real motion of the iceberg with the associated error matrixs. During the forecasting, the filter runs without inputs, ie without measured positions. Chapter 4. Application of the Kalman filter to iceberg motion 30 Figure 4.4: Real iceberg trajectories Chapter 4. Application of the Kalman filter to iceberg motion 31 Let xp be the last filtered position. Without measured position at time tp+i. xp+i = xp+1 = Axp. More generally, Xp+n = A" xp. The process equation gives: xp+i = Axp + wp At the following step, xp+2 = A(Axp + wp) + wp+i and so forth. n Therefore, xp+n = Anxp + ^ A^Wp^ (4.32) «=i = An is the process matrix associated with the n-step forecasting, and u,(n) = ^An-,Wp+._i (4 33) i=l is the associated noise matrix. The error matrix associated with the n-step focasting is Pp-+n = E[{xp+n - x;+n)(xp+n - x-+n)T] (4.34) Having forecasted positions and errors, it is now possible to evaluate the risk of impact. Chapter 5 Risk of impact calculation 5.1 A numerical algorithm The following risk calculation procedure has been developed by A.B. Dunwoody [11]. It is not the purpose of this papre to explain the complete scheme, but the general steps should be presented. Let n(T) be the number of impacts between time=0 and T. Supposing that p(n(T) > 0) <C 1, n(T) can be integrated over the probability of impact between small intervals (t,t -f St) as follows: p(n(tf,t + 8t) = 1) is integrated over all possible velocities upon impact. Let u(t) be the velocity of the iceberg at the time of impact. Then: (5.35) p(n(t, t + 6t) = l)= / p(u(t) = Ukn(t, t + St) = l)dU (5.36) o The joint probability of impact and speed is: 32 Chapter 5. Risk of impact calculation 33 p(u(t) = Ukn(t, t + 6t) = l) = rV* ra+-Jo Ja+% ±U2cos(d-a)6tx ( u(t) = U cos 0 sin0 cos a sin a I > d9da (5.37) with d being the distance between iceberg and platform upon impact and r(t) = (x(t), y(t)) being the position of the ice feature at time t. The joint probability distribution of the position and velocity is a multivariate normal density function, because the physical system is linear and excited by gaussian white noise. p(r(t) = Rku(t) = U) = (27r)2|D|5 6XP^ 2 X - E[x] " ( , X - E[x] ^ X -E[x] zr1 X -E[x] )(5.38) y -E[y] y ~E[y] \y - m, { V - E[y] j with D = Di (0) \(0) D2 E[(x - E[x])2] E[(x - E[x])(x - E[x])] \ E[(i - E[x]){x - E[x]) E[(x - E[x])2) ( = / E[(y - E[y})2} E[(y - E[y])(y - E[y))} { E[(y - E[y})(y - E[y\) E[(y - E[y})2} Chapter 5. Risk of impact calculation 34 where (t) has been omitted to save writing. is indeed the forecasting of the process vector at time t, 1 E[x) ^ E[x] E[y] \ E[y] j and D is the mean square error matrix previously noted as P~. Risk calculation is straightforward using the forecasted data provided by the Kalman filter. In risk integration, the time increment is the time separating two forecasted po sitions. At each step into the future, the joint probability distribution of position and velocity is the normal function of process vector and error matrix at the corresponding time. During the forecasting, the Kalman filter gives positions and error matrices. The risk calculation procedure can be implemented. The risk calculation includes three levels of integration: Level 1 Integration over time of the risk of impact (equation 5.36). Level 2 At each time step, integration over all possible drift speeds of the joint proba bility of a particular speed and impact (equation 5.37). Level 3 Calculation of the joint probabihty of speed and impact as an integral of the joint distribution of position and velocity over an area surrounding the platform (equation 5.38). The three levels are calculated using Simpson's rule. Because of the possible accu mulation of errors, it is necessary to be careful with integration domains. Uncertainties at level 3 are highly amplified through the second and first levels. Chapter 5. Risk of impact calculation 35 5.2 Application of the filter to a real trajectory The following example will describe step by step the application of the Kalman filter to the motion of an iceberg. Figure (5.5) shows an iceberg trajectory observed on the east coast of Canada. Po sitions were plotted relative to the platform at hourly intervals. The first study of this trajectory consists of an estimation of the process noise applied on velocity. Velocities and accelerations are deduced from positions. The process noise is estimated as the stan dard deviation of acceleration. In this case, = 0.15 and -Efu^] = 0.31 (units are km2h~4). The measurement uncertainty is taken as 400 meters, ie E[v2] = E[v2] — 0.16 (units are km2). The treatment of the trajectory starts at time = 2 hours, where the process vector is taken as the vector of observed position and velocity. During the treatment of the trajectory, the procedure is repeated as follows: for instance, let JC4 be the last filtered position at time 4. The forecast of the process vector at time 5 is x$ = Ax A- This is the one step predicted position found in figure (5.6). The observed position at time 5 is then obtained from the radar z5 (figure (5.7)) and the last filtered position is i5 = xl + K5(z5 — C5xj) (figure (5.8)). It requires 5 steps for gain and errors to converge towards stationary values. Figure (5.9) shows the curves of observed, step by step predicted, and filtered trajec tories. The filter tends to smoothen the sudden changes in the iceberg's velocity. If the noises are taken smaller, the gain of the filter is higher, and therefore the weight given to observation when calculating the filtered value is higher, too. Following the last observed position, the filter is running without input. The filter forecasts the future trajectory (figure (5.10)). The Root Mean Square (r.m.s) errors on position can be used as a measure of the Chapter 5. Risk of impact calculation 36 -7 -r y -12 -11 -10 -9 -8 -7 -6 x axis (km) Figure 5.5: Actual iceberg trajectory used in example Chapter 5. Risk of impact calculation 37 • Observed trajectory • Filtered trajectory £ One step forecasting 1 1 fi 1 ° -8 -7 forecasting for time=5 \_. || first forecasting g 2 \ first filtered position^ jj -9 -• time=0 first observed position -10 -5 Figure 5.6: One-step-ahead forecasted position Chapter 5. Risk of impact calculation 38 Figure 5.7: The newly observed position Chapter 5. Risk of impact calculation 39 -8.5 -7.5 -6.5 last filtered position ii observation at time=5 a x i s • -7 — -8 -10 x axis (km) • Observed trajectory aV One step forecasting • Filtered trajectory Figure 5.8: Filtered position Chapter 5. Risk of impact calculation 40 Figure 5.9: Filtered trajectory Chapter 5. Risk of impact calculation 41 Figure 5.10: Forecasted trajectory Chapter 5. Risk of impact calculation 42 uncertainty over position in the forecasting. It is the square root of the variance of the error: / E((x(t) - x(t))2) \ { E((y(t) ~ Vi*))3) I Ellipses can be drawn, with the two r.m.s position errors as radii. Supposing that the errors are gaussian, there is approximatively a 42% chance that the iceberg will be within the ellipse at the corresponding time. Errors on positions start to grow after the last known position (figure (5.11)), and become rather important after 20 hours of prediction (figure (5.12)). 5.3 An example of risk calculation The idea is to compare the results of the risk calculation, and those given by a simulation. An iceberg trajectory has been created to serve as the starting point of the test. The noise variances -E[tuJ] and .Efu^] are given and are both equal to 0.1 km2h~4. The trajectory is developed from a given initial position and velocity, using the process equation, xk+1 - Axk + w w is randomly chosen at each time step within a gaussian distribution with zero mean and £[to2] variance. Because the trajectory only has 26 points, the variances are not strictly equal to 0.1 km2h,-4, but E[w2) = 0.096 km2h-4 and E[wl) = 0.081 km2h~4. This is the real trajectory shown in figure (5.13). In order to create conditions as close as possible to reality, the trajectory should be corrupted with measurement noise. The measurement noise is taken as 250 meters. The "observed" trajectory is created by randomly taking each of its points within a gaussian Chapter 5. Risk of impact calculation 43 Figure 5.11: Increase of uncertainty perimeters (r.m.s errors) Chapter 5. Risk of impact calculation 44 Figure 5.12: Increase of position uncertainty with time Chapter 5. Risk of impact calculation 45 Figure 5.13: A computed trajectory Chapter 5. Risk of impact calculation 46 distribution (figure (5.14)). The variance of the distribution is the given measurement noise. In order to implement the Kalman filter, an estimate of the process noise is required from the observed trajectory. Velocities and acceleration of the ice feature are deduced from the positions. The variances of acceleration are then calculated and given i^[^i] = 0.11km2h~4 and E[wl\ = 0.09 km2h~4. This is the value of the noise given to the filter, the only information available to the forecaster. In this case, and because the measurement noise is smaller than the process noise, the measurement of the process noise over the observed trajectory gives a good level of accuracy. If the filter is run several times, using at each loop a new noise calculated with the last filtered trajectory, we obtain after convergence 0.024 and 0.019km2h~A. The filter tends to make the trajectory "smoother" and reduces the noise to a value close to the measurement noise. Therefore, when the process noise is higher than the measurement noise, which is true in most of the cases, the first estimation of the process noise is better than a value calculated with the filter. The final filtered trajectory, as well as 25 hours of forecasted positions and the r.m.s error circles are shown in figure (5.15). The platform is situated at the point (0,0) of the figure. The impact perimeter is taken as 1 km, ie there would be impact if the center of the iceberg would come closer than 500 meters to the center of the platform. According to the forecasted path of the iceberg, the platform is most vulnerable to impact from the iceberg 14 hours after the last observation (figure (5.16)). The risk calculations involve certain precautions which merit discussion, if realistic results are to be obtained. The time steps at level 1 are given by the time lag between two forecasted positions. In order to obtain a good numerical integration, the r.m.s error circles must overlap, or in other words, the time lag between positions should be small enough to observe a Chapter 5. Risk of impact calculation 47 Figure 5.14: Corresponding "observed" trajectory Figure 5.15: Trajectory forecasting and position error perimeters Chapter 5. Risk of impact calculation Figure 5.16: Forecasted positions and corresponding times Chapter 5. Risk of impact calculation 50 progressive change in the risk of impact with time. If it is not, a solution could be to interpolate some values of risk, or reduce the time step in the forecasting algorithm. At level 2, the velocity range has to be defined for integration. A possible solution is to calculate the average velocity of the iceberg during filtering. The hmit of the velocity domain is then taken as n times this velocity. A longer but more accurate solution is to numerically "visualize" the joint probability of velocity and impact, and then to define at each time step a velocity maximum for integration. The angles a and 9 appearing in equation 5.36 have to be chosen small enough. Tests have shown a convergence of the risk integration in most of the cases at a — 0 — 10 deg. In the presented case, the complete numerical integration gives a risk of 7% . Figure (5.17) shows the distribution of risk of impact versus time. The highest risk is obtained at time = 10 hours, 4 hours before the expected time of impact. This difference is due to the increase of uncertainties in positions and velocities with time. At the expected time of impact, the r.m.s position error radius is equal to 10 km. If the errors are supposed to be gaussian, the iceberg has about a 42 percent chance of being in this circle. At time =11 hours, the radius of the circle is 7 km. A complete discussion concerning risk calculation can be found in [11]. 5.4 Comparison of risk calculation and simulation The principle of the simulation is to run random trajectories from the last real position. The risk is the number of recorded impacts divided by the number of trajectories. The random trajectories are created using the process equation and the real process noise, both components being equal to 0.1 km2h~4. In this case, 416 impacts have been recorded for 5000 trajectories, giving a risk of 8.4%. The distribution of time of impact is given in figure (5.18). Although the shape of Chapter 5. Risk of impact calculation 51 0.012 T time of forecasting (hours) Figure 5.17: Calculated risk distribution Chapter 5. Risk of impact calculation 52 the curve is similar to the one calculated in the previous section, the risk is 20% higher. A risk difference of 10-30% has been recorded between calculations and simulations, for different initial "real" trajectories. The reason for this difference is that the sim ulation doesn't take into account the uncertainties over the last known position. This position is the starting point of the random trajectories. The simulation only considers the dispersion of the trajectories due to the process noise. There are already uncertainties concerning position and velocity of the last filtered value. The situation at this particular point of the trajectory requires careful analysis. In the developed example, real and filtered process values are very close to each other. While the position difference lies at about 5%, the velocity difference equals approximately 15%. The r.m.s position error is 0.184 km, the r.m.s error for velocity is 0.37 km . Practically, this means that the direction of the velocity vector could vary by ±35deg relative to the filtered value, within a 90% confidence limit. Figure (5.19) shows the corresponding dispersion of straight trajectories starting from this point. In order to compare results of risk calculation and simulation, equivalent conditions should be considered in the simulation. The importance of position uncertainty is small. The velocity uncertainty can be considered by running several sets of random trajectories with different initial velocities. The initial range of velocity has been divided into 13 values on both x and y axes. The risk is recorded for each of the 169 tests. The final risk is the integral over velocities on x and y of the corresponding risks. This operation is highly time consuming. If every experimental set contains 200 random trajectories, 33,800 trajectories have been created in the end. The resulting distribution of time of impact is given in figure (5.20). The overall impact risk is 6.8%, which should be compared with the 7% result of the risk calculation algorithm. Figures (5.17) and (5.20) are also very similar. Chapter 5. Risk of impact calculation 53 0.018 T 0.016 p 0.014 --r 0 0.012 + b a 0.01 b 1 0.008 4-1 i 0.006 + t y 0.004 --0.002 --0 9 11 13 15 17 19 21 23 25 forecasting time (hours) Figure 5.18: Risk distribution obtained with a simple simulation! Chapter 5. Risk of impact calculation 54 Figure 5.19: Uncertainty over the direction of the future trajectory Chapter 5. Risk of impact calculation 55 0.012 x 0.01 --P forecasting time (hours) Figure 5.20: Risk distribution obtained with a complete simulation! Chapter 6 Conclusion The Kalman filter is a powerful tool for treating noiselike signals. Iceberg motions are particularly noisy: measured trajectories show high unpredictability. Furthermore, this application of the Kalman filter can easily be used with the risk calculation algorithm developed by B.Dunwoody [11]. The practical advantages are multiple: 1. The most important feature is certainly the hmited amount of information required to run the prediction algorithm, and therefore the risk calculation. The user needs only the measurement uncertainty of the radar and the observed trajectory of the iceberg in order to estimate the risk of impact to the platform. 2. The Kalman filter is recursive. Its parameters are self-adaptive to new inputs. The filter keeps track of the past. A new value of the risk can be easily calculated as soon as a new measurement of the position of the iceberg is available. Only the process noise has to be recalculated after a certain period of time in order to reach a good level of accuracy. Several trajectories can be treated with Hmited computations. 3. The filter does not have to be modified to accept deterministic input. Its param eters, such as gain or errors, are not changed by deterministic input. They can be the influence of a mean current or wind. A deterministic equation according to the known influence can simply be added to the filter. It could improve the predicted trajectory without changing errors over expected position and velocities. This might be of major 56 Chapter 6. Conclusion 57 interest in the vicinity of the platform, or in areas with known steady flow or wind fields. 4. The time step between observations can be modified. Observations are usually made more frequently as the iceberg approaches the platform. The Kalman filter is a flexible and adaptive method serving as a new tool in the field of iceberg motion forecasting. Some improvements of the actual version have to be considered: 1. Time scale of observations: the closer in time the observations are, the better the process noise can be calculated, and the more errors due to measurement uncer tainties can be reduced. For instance, two opposite sided errors on the last two positions can create a highly biased direction for the forecasted trajectory. More frequent observations will reduce this risk. 2. The filter gives the "best" estimate of the process vector according to the given parameters. The system transition matrix, which links two steps of the process, should be numerically improved, in order to give a better image of the reality. This could involve discrete time modeling of the iceberg motions mechanism, as well as the study of recorded trajectories for a given site. Bibliography [1] A. V. Balakrishnan, "Kalman filtering theory" Optimization Software, Inc. Publi cation Division, New York, 1984. [2] P. Ball, H. S. Gaskill and R. J. Lopez, "An analysis of two data sets collected at drill sites in the Labrador sea" C-Core Publication, Memorial University of Newfoundland pp. 81-82, 1981. [3] P. Ball, H. S. Gaskill and R. J. Lopez, "Environmental data requirements for a real time iceberg motion model"pp. 1369-1378, 1981. [4] S. M. Bozic, "An introduction to digital and Kalman filtering" Edward Arnold, 1979. [5] R. G. Brown, "Introduction to random signal analysis and Kalman filtering" John Wiley & Sons, editors, 1983. [6] D. S. Davison and R. McKenna, "Statistical modelling techniques for iceberg mo tion in West Baffin Bay" Atmosphere-Ocean, 16th Annual Congress, University of Ottawa, May 1982. [7] R. T: Dempster, "The measurement and modeling of iceberg drift" IEEE Ocean '74 pp. 125-129,1974. [8] R. T. Dempster, "Characteristics of iceberg mechanics" Proceedings Symp. on the Physics and Mechanics of Ice, International Union of Theoretical and Applied Me chanics, Copenhagen, Denmark pp. 38-50, 1979. 58 Bibliography 59 [9] R. T. Dempster, "Tacticaliceberg management", preprints of a Seminar on Iceberg Management in Offshore Exploration, Production and Transportation, Memorial University of Newfoundland, Extension Services, Continuing Engineering Educa tion, St. John's, 1982. [10] R. T. Dempster and D. S. Sodhi, "Motion of icebergs due to changes in water currents" IEEE Ocean '75 pp. 348-350, 1974. [11] A. B. Dunwoody, "Operational risk evaluation of an ice feature impact" Journal of Offshore Mechanics and Arctic Engineering, pp. 96-101, February 1990. [12] C. Garrett, "Statistical prediction of iceberg trajectories" Cold Regions Science and Technology pp. 255-266,1985. [13] C. Garrett, M. Hazen, F. Majaess and J. F. Middleton, "Analysis and prediction of iceberg trajectories" Dalhousie University, Halifax, Feb. 1985. [14] H. S. Gaskill and J. Rochester, "A new technique for iceberg drift prediction" Cold Regions Science and Technology pp. 223-234,1984. [15] H. S. Gaskill and B. Terry, "Tactical iceberg management - drift forecasting", preprints of a Seminar on Iceberg Management in Offshore Exploration, Production and Transportation, Memorial University of Newfoundland, Extension Services, Continuing Engineering Education, St. John's, 1982; [16] R. E. Kalman, "A new approach to hnear filtering and prediction problems" ASME-Joumal of Basic Engr. pp. 35-45, 1960. [17] B. W. Leach, "An introduction to Kalman filtering" National Research Council Canada, 1984. Bibliography 60 [18] D. G. Mountain, "On predicting iceberg drift" Cold Regions Science and Technology pp. 273-282,1980. [19] D. B. Muggeridge, N. P. Riggs and K. Shirasawa, "The drift of a number of idealized model icebergs" Cold Regions Science and Technology pp. 19-30, 1984. [20] J. E. Murray, "The drift, deterioration and distribution of icebergs in the North Atlantic ocean" Ice Seminar, Canadian Institute of Mining and Metallurgy, 1969. [21] S. D. Smith and E. G. Banke, "A numerical model of iceberg drift" Proceedings of the 6th International Conference on Port and Ocean, Engineering Under Arctic Conditions, Quebec, 1981. [22] S. D. Smith and E. G. Banke, "The influence of winds, currents and towing forces on the drift of icebergs" Cold Regions Science and Technology pp. 241-255, 1983. [23] M. El-Tahan and D. S. Sodi, "Prediction of an iceberg drift trajectory during a storm" Annals of Glaciology pp. 77-82, 1980. [24] M. El-Tahan, "Iceberg trajectories, C-CORE", personal communication, 1990. [25] N. Wiener, "Extrapolation, interpolation, and smoothing of stationary time series" Wiley, 1949. Appendix A : Flowchart of the filtering software Pinit, Q, C, R, A 1' Initial erroi Pi covariance nit / Predicted covariance Pk~ = A^P^A^ + Qk-y Initial estimate of the process vector Prediction *fc = Ak-iXk-i Filter gain Kk = Pk-Cl{CkPk~Cl + Rk)-Estimate of present output ifc = Ckxk Error covariance for updated estimate Pk = (I- KkCk)Pk~ Error efc = zk - *k Correction Kkek Current estimate xk = xk + Kkek \j, ^^rerun* -<-<^for another^ Present output /ru Store last error covariance Piatt Store last estimate of the process vector Xlatt 61 Appendix B : Flowchart of the forecasting software Forecasting n steps ahead Last error covariance Piatt Last estimate of the process vector Process noise covariance calculation Prediction xn = AnX(att Predicted covariance P~ = AnPla,t(An)T + Storage of P and x 62 Appendix C : Flowchart of the risk calculation software Joint probability distribution of position and velocity p(r(t) = Rku{t) = U) = 1 (2*)2|P- fyexp("2 fx — E[x] \ T ( x - E\x) \ x - E[x] x - E[x] V - E[y] y - Ely) \y- E[y] ) \ i/ - E[y] J Joint probability of impact and speed 2» ,0+^ p(u{t) = Ukn{t,t + St) = 1) = - / * /° : JO ' fl72cos(0- a)Stx u(t) = u{ COsa ) dOda Probability of impact during small time interval (between two steps) p(n(t, t + 6t) = l)= (°° p(u{t) = Ubn(t, t + 6t) = l)dU Jo * Probability of impact 63
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of the Kalman filter to iceberg motion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of the Kalman filter to iceberg motion forecating Simon, Christophe 1990-09-27
pdf
Page Metadata
Item Metadata
Title | Application of the Kalman filter to iceberg motion forecating |
Creator |
Simon, Christophe |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | The objective of this study is to develop an application of the Kalman filter for filtering and forecasting iceberg positions and velocities in order to calculate the risk of impact against a fixed structure or stationary vessel. Existing physical and probabilistic models are reviewed. Physical models are essentially based on the response of the iceberg to the forces acting on it. Statistical models forecast the motion of the iceberg based on past observations of the trajectory. A probabilistic iceberg trajectory model is used in this study so that uncertainties in the trajectory forecast can be explicitly included. The technique of Kalman filtering is described and applied to forecast future positions and velocities of an ice feature, including uncertainties.The trajectory forecast combined with a risk calculation, yields the probability of impact and the probability distribution of the time until impact. Numerical simulations are provided to demonstrate and validate the procedure. |
Subject |
Kalman filtering Icebergs -- Forecasting |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0098011 |
URI | http://hdl.handle.net/2429/28723 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1990_A7 S55.pdf [ 2.94MB ]
- Metadata
- JSON: 831-1.0098011.json
- JSON-LD: 831-1.0098011-ld.json
- RDF/XML (Pretty): 831-1.0098011-rdf.xml
- RDF/JSON: 831-1.0098011-rdf.json
- Turtle: 831-1.0098011-turtle.txt
- N-Triples: 831-1.0098011-rdf-ntriples.txt
- Original Record: 831-1.0098011-source.json
- Full Text
- 831-1.0098011-fulltext.txt
- Citation
- 831-1.0098011.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0098011/manifest