Discrete-time Kalman filter for heave motion estimation

Joel Reis, Pedro Batista, Paulo Oliveira, and Carlos Silvestre

Ocean Engineering, vol. 277, 114240, June 2023

https://doi.org/10.1016/j.oceaneng.2023.114240

Accepted Version
Level of access, as per info available on SHERPA/ROMEO
http://www.sherpa.ac.uk/romeo/search.php
Discrete-time Kalman filter for heave motion estimation⋆
Joel Reisa , Pedro Batistab , Paulo Oliveirab,c and Carlos Silvestrea,b,∗,1
a Faculty of Science and Technology, University of Macau, Macau, China
b Institutefor Systems and Robotics, Instituto Superior Técnico, Universidade de Lisboa, Lisboa, 1049-001, Portugal
c Associated Laboratory for Energy, Transports and Aeronautics, Institute of Mechanical Engineering, Instituto Superior Técnico, Universidade de

Lisboa, Lisboa, 1049-001, Portugal

ARTICLE INFO Abstract

Keywords: This paper addresses the problem of estimating the heave motion of a platform using biased
Heave motion estimation measurements of an accelerometer. We develop a general framework wherein convenient and well-
Kalman filter known properties of trigonometric functions are exploited to devise a linear system whose state
Observability encompasses implicit representations of wave amplitudes and phase shifts, as well as a constant
Sea waves spectrum sensor bias. The observability of the system is analyzed, followed naturally by the implementation
of a discrete-time linear time-varying Kalman filter with global asymptotic stability guarantees.
Our proposed methodology is validated with realistic numerical examples, including an accurate
representation of a continuous wave spectrum within an ocean context.

1. Introduction to heave estimation concern measurements of significant

wave height and sea surface elevation, see, e.g., Collins,
The offshore industry is the scene of complex engineer- Lund, Waseda and Graber (2014), Nielsen, Brodtkorb and
ing operations with an intricate relationship between the
Sørensen (2019) and references therein.
surface and the underwater and aerial media. For instance,
Attaining a high degree of accuracy in heave informa-
towards the elimination of safety and operational hazards, tion is possible by using real-time kinematic positioning
predicting and monitoring the heave of a surface vessel is to correct for errors in global navigation satellite systems
critical during: i) the installation of structures and facilities (GNSSs), see, e.g., Blake (2007) and Kielland and Hagglund
in a marine environment, e.g., bridges and underwater labs; (2015). However, differential GNSS may not always be
ii) the remote control of subsea vehicles connected via an
available and, in a worst case scenario, even GNSS alone
umbilical cord; and iii) the vertical take-off and/or landing of
could perform erratically. The reliance on other sensors
aircraft. The available literature on this urgent topic includes onboard the surface vessel proves not only practical but also
a large and heterogeneous collection of articles; in addition advantageous in terms of operational redundancy.
to the relevant references cited in this paper, the interested Using compact inertial measurement units (IMUs) to
reader is also referred to, for example, Johansen, Fossen, estimate the heave turns out to be rather cheap and efficient,
Sagatun and Nielsen (2003), Skaare and Egeland (2006),
as demonstrated in the works of Küchler, Eberharter, Langer,
Woodacre, Bauer and Irani (2015) Li, Wei, Guo and Zhu
Schneider and Sawodny (2011a), Küchler, Mahl, Neupert,
(2017), and references therein, for a more general perspec- Schneider and Sawodny (2011b), and Auestad, Gravdahl and
tive of the problem. In particular, because sea waves cause Fossen (2013). We note that a similar strategy is pursued by
large relative displacements of deck cranes, much attention the authors in this article but without resorting to lineariza-
and effort have been put into the industrial development tions or implicit observations of the heave. The work in Li,
of active heave compensation winch systems that can act
Li and Wen (2022) is particularly noteworthy for its solution
swiftly in the presence of mildly rough seas (see, e.g., Mac-
based on two IMUs, and the work in Richter, Schaut, Walser,
Gregor (2022), Wago (2022)). Other applications related Schneider and Sawodny (2017) for its comprehensive active
⋆
This work was supported in part by the Macau Science and Technol- heave compensation approach that includes estimation, pre-
ogy Development Fund under Grant FDCT/0146/2019/A3, in part by the diction and control methods, with the former obtained by
University of Macau under Grant MYRG2020-00188-FST, and in part by applying adaptive filtering methods.
the Fundação Para a Ciência e a Tecnologia (FCT) through LARSyS–FCT
Plurianual funding 2020-2023 and LAETA-FCT Project under Grant A method based on a coarse- and fine-tuning fixed-grid
UIDB/50022/2020 and through the FCT project DECENTER under Grant wavelet networks is presented in Huang, Jiang and Zou
LISBOA-01-0145-FEDER-029605, funded by the Programa Operacional (2021) for online prediction of the coupled heave-pitch mo-
Regional de Lisboa 2020 and PIDDAC Programs. The work of Joel Reis tions of a ship in irregular waves. The prediction of the heave
was supported by the FDCT Funding Scheme for Postdoctoral Researchers
of Higher Education Institutions 2021 under Grant 0030/2021/APD. motion of a vessel at different sea states is also studied by
∗ Corresponding author Tang, Yao, Li, Wang and Hu (2022), where, inspired by the
joelreis@um.edu.mo (J. Reis); pbatista@isr.tecnico.ulisboa.pt (P. flexible and powerful optimization capability of the particle
Batista); pjcro@isr.tecnico.ulisboa.pt (P. Oliveira); csilvestre@um.edu.mo swarm optimization algorithm, the authors propose a hybrid
(C. Silvestre)
ORCID (s): 0000-0002-9281-6829 (J. Reis); 0000-0001-6079-0436 (P.
autoregressive moving average model. Also aiming at en-
Batista); 0000-0002-5799-390X (P. Oliveira); 0000-0002-5096-5527 (C. suring high accuracy results under unknown sea conditions,
Silvestre) the work in Ben, Gao, Wei, Gong and Li (2022) employs
1 On leave from Instituto Superior Técnico, Universidade de Lisboa.

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 1 of 9

 Discrete-time Kalman lter for heave motion estimation
a band-pass filter instead of the high-pass filter existing in (or vector) of zeros; if the matrix is square of dimension
traditional methods; more specifically, the authors develop 𝑚 × 𝑚, then, for notational simplicity, we just write 𝟎𝑚 .
an adaptive frequency estimator based on the weighted- Similarly, the symbol 𝐈𝑚 denotes an 𝑚 × 𝑚 identity matrix.
frequency Fourier linear combiner method. The transpose operator is denoted by the superscript (∙)𝖳 . A
For the case when observations of the environment in- positive definite matrix 𝐌 is denoted by 𝐌 ≻ 𝟎. The special
clude the position of a platform and the surface elevations orthogonal group of order three is denoted by 𝖲𝖮(3) ∶=
obtained from wave gauges, Xie, Zhao and Luo (2022) {𝐗 ∈ ℝ3×3 ∶ 𝐗𝐗𝖳 = 𝐗𝖳 𝐗 = 𝐈3 , det(𝐗) = 1}. For conve-
propose an extended deep reinforcement learning framework nience, 𝐞{𝑛},𝑗 ∈ ℝ𝑛 is an 𝑛-dimensional unit column vector
to find the optimal control strategy to dissipate wave energies with all elements zero except for the element at index 𝑗 ≤ 𝑛,
with an active-controlled heaving horizontal plate. In harsh which is one, e.g., if 𝑛 = 5, then 𝐞{5},4 = [0 0 0 1 0]𝖳 . In order
sea conditions, the safety and efficiency of subsea operations to simplify expressions, whenever the argument of a scalar
can be maximized by regulating the length of an umbilical operator is a vector, said operator is applied, in an element-
cable through an active heave compensation system. Zhou, wise fashion, to each entry
( )of its(vector
) argument,
( ) e.g., given
Cao, Yao and Lian (2021) develop a robust hierarchical 𝐱 ∈ ℝ𝑛 , sin(𝐱) ∶= [sin 𝑥1 sin 𝑥2 … sin 𝑥𝑛 ]𝖳 ∈ ℝ𝑛 .
multi-loop control scheme that combines nonlinear model
predictive control with integral sliding mode control; the
heave motion is obtained by data fusion of the heave mo- 3. Problem statement
tion measured through an motion reference unit, which can Let the vertical position of a heaving platform be de-
directly feedback the vessels’ motion state through a twice scribed by the following time-dependent finite sequence:
integration of the measured vector acceleration. Naturally,
smoothing out data should always be taken into considera- ∑
𝑁
( )
𝑧(𝑡) ∶= 𝐴𝑖 sin 𝜔𝑖 𝑡 + 𝜙𝑖 ∈ ℝ, (1)
tion, as convincingly argued in El-Hawary (1982), where a
𝑖=1
model for the heave dynamics is postulated that provides the
basis for formulating the heave extraction problem as one of where 𝑁 ≥ 1 represents the number of principal waves,
Kalman filtering in the theory of optimal linear estimation, and where 𝐴𝑖 , 𝜔𝑖 , 𝜙𝑖 ∈ ℝ correspond to the amplitude,
with the filter ultimately reducing the residual heave effects. frequency and phase, respectively, of each principal wave. In
As opposed to Küchler et al. (2011a), where an extended this work, it is assumed that all three wave parameters, i.e.,
Kalman filter is implemented to estimate the heave motion 𝐴𝑖 , 𝜔𝑖 , and 𝜙𝑖 are constant. Furthermore, it is assumed that:
with high accuracy, and as opposed to Auestad et al. (2013), all 𝑁 frequencies are positive and different, i.e., 𝜔𝑖 ≠ 𝜔𝑗 > 0
where a continuous model is discretized using the one step for 𝑖, 𝑗 ∈ {1, 2, … , 𝑁}, 𝑖 ≠ 𝑗; amplitudes are positive;
Euler method, in this work we propose, inspired by the phases are set in the interval [−𝜋, 𝜋]. The derivative of (1),
results reported in Batista, Oliveira and Silvestre (2018), a which represents the heave speed, is given by
solution that exploits the sine addition formula to devise
a discrete-time linear time-varying (DT-LTV) system that ∑
𝑁
( )
̇ =
𝑧(𝑡) 𝜔𝑖 𝐴𝑖 cos 𝜔𝑖 𝑡 + 𝜙𝑖 ∈ ℝ. (2)
avoids linearized approximations, offering guarantees of
𝑖=1
global asymptotic stability. It is said formula that allows us
to set this work apart from Küchler et al. (2011b), where Consequently, the coordinated acceleration associated with
a linear discrete-time, but time-invariant Kalman filter is the oscillatory behavior defined in (1) is calculated as
employed to estimate the vertical motion of a crane tip
and adapt the parameter vectors by comparing the mod- ∑
𝑁
( )
̈ =
𝑧(𝑡) 𝛽𝑖 sin 𝜔𝑖 𝑡 + 𝜙𝑖 ∈ ℝ, (3)
eled motion with the measured motion. Similarly, a linear
𝑖=1
discrete-time Kalman filter is also employed in Tang et al.
(2022), but it is solely used to filter the collected acceleration with constant 𝛽𝑖 ∶= −𝐴𝑖 𝜔2𝑖 < 0, for 𝑖 = 1, 2, … , 𝑁.
signals. Moreover, we show that our strategy can also work Suppose now that the specific acceleration of the plat-
in relatively unknown sea conditions, granted that a region form is periodically sampled with sampling frequency 𝑓𝑠 ∶=
of the wave spectrum, which is realistically and practically 1∕𝑇𝑠 ∈ ℝ, where 𝑇𝑠 > 0 is the sampling time. These
attainable. Like El-Hawary (1982), our Kalman filter may samples are assumed to be collected from a bias-corrupted
be applied to post-experiment data records, or in a real-time accelerometer, whose readings we express here as
operation working directly on measurements collected from
an IMU. Moreover, our realistic simulations featuring sensor ̈ 𝑘 ) + 𝑔𝐞𝖳{3},3 𝐑(𝑡𝑘 )𝐞{3},3 + 𝑏 ∈ ℝ,
𝑎𝑚 (𝑡𝑘 ) ∶= 𝑧(𝑡 (4)
noise and an offset allow for a exhaustive communication of
the validation and assessment results. where 𝑡𝑘+1 ∶= 𝑡𝑘 + 𝑇𝑠 represents the sampling instant,
𝐑 ∈ 𝖲𝖮(3) denotes the attitude of the platform, 𝑔 ∈ ℝ is
the local constant acceleration of gravity, and 𝑏 ∈ ℝ is an
2. Notation unknown accelerometer bias, also assumed constant.
Throughout the paper, a bold symbol stands for a multi- The main objective of this work is to obtain an estimate
dimensional variable. The 𝑛-dimensional Euclidean space is of 𝑧(𝑡𝑘 ) using biased accelerometer measurements as de-
represented by ℝ𝑛 . The symbol 𝟎𝑚×𝑛 denotes an 𝑚×𝑛 matrix scribed by (4). More specifically, and based on Kalman filter

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 2 of 9

 Discrete-time Kalman lter for heave motion estimation

Table 1 We then create an auxiliary column vector containing all

List of main symbols used throughout this paper known frequencies, given by
Symbol Description [ ]𝖳
𝝎 ∶= 𝜔1 𝜔2 … 𝜔𝑁 ∈ ℝ𝑁 . (11)
𝑁 number of principal waves, 𝑖 = 1, 2, … , 𝑁
𝜔𝑖 angular frequency of 𝑖-th wave The complete system state vector may thus be defined as
𝜙𝑖 phase shift of 𝑖-th wave
𝑥1,𝑖 invariant auxiliary system state of 𝑖-th wave [
𝑥2,𝑖 invariant auxiliary system state of 𝑖-th wave 𝐱[𝑘] ∶= 𝑥1,1 [𝑘] 𝑥1,2 [𝑘] … 𝑥1,𝑁 [𝑘]
𝑏 constant accelerometer bias ]𝖳 (12)
𝑔 constant gravitational acceleration 𝑥2,1 [𝑘] 𝑥2,2 [𝑘] … 𝑥2,𝑁 [𝑘] 𝑏[𝑘] ∈ ℝ2𝑁+1 ,
𝑇𝑠 sensor sampling time
for which we can express a DT-LTV system given by
{
theory, we set out to obtain estimates of the sensor offset 𝑏, 𝐱[𝑘 + 1] = 𝐀𝐱[𝑘], (13a)
followed by an implicit determination of all wave amplitudes 𝑦[𝑘 + 1] = 𝐂[𝑘 + 1]𝐱[𝑘 + 1], (13b)
and phase shifts that enables the reconstruction of the heave
motion characterized by (1). where the time-invariant matrix 𝐀 ∈ ℝ(2𝑁+1)×(2𝑁+1) is
defined as an identity matrix given by
4. DT-LTV system design ⎡ 𝐈𝑁 𝟎𝑁 𝟎𝑁×1 ⎤
The (nominal) derivations presented in this section are 𝐀 ∶= ⎢ 𝟎 𝐈𝑁 𝟎𝑁×1 ⎥ ≡ 𝐈2𝑁+1 , (14)
⎢ ⎥
based on the fact that the wave frequencies are known. ⎣𝟎1×𝑁 𝟎1×𝑁 1 ⎦
Reasonably, such knowledge may be obtained a priori via,
e.g., a fast Fourier transform of the sampled signal (4). Also
noteworthy is that knowledge of phase or amplitude is not 𝑦[𝑘] ∶= 𝑎𝑚 [𝑘] − 𝑢[𝑘] ∈ ℝ, (15)
required, which represents a much less demanding aspect
from an implementation point of view in contrast with other and, finally, where, according to (10), (11) and (15), 𝐂[𝑘] ∈
methods in the literature. ℝ1×(2𝑁+1) is constructed as
Table 1 contains a list of the main symbols used through- [ ( ) ( ) ]
out the paper. 𝐂[𝑘] ∶= sin 𝝎𝑘𝑇𝑠 𝖳 cos 𝝎𝑘𝑇𝑠 𝖳 1 . (16)
We proceed by noticing that, in view of the angle ad-
dition theorem of the sin(∙) function, we can write, for 𝑖 = Assumption 1. The sampling frequency 𝑓𝑠 verifies the
1, 2, … , 𝑁, the following identity: Nyquist sampling rate for the largest wave frequency. This
( ) ( ) ( ) means that, without loss of generality, given 𝜔𝑁 > 𝜔𝑖 , for
𝛽𝑖 sin 𝜔𝑖 𝑡𝑘 + 𝜙𝑖 = sin 𝜔𝑖 𝑡𝑘 𝛽𝑖 cos 𝜙𝑖 𝑖 = 1, 2, … , 𝑁 − 1, we have
( ) ( ) (5)
+ cos 𝜔𝑖 𝑡𝑘 𝛽𝑖 sin 𝜙𝑖 . 𝜔𝑁
𝑓𝑠 > 2𝑓𝑁 , with 𝑓𝑁 ∶= . (17)
For notational simplicity, we define 2𝜋

𝑢[𝑘] ∶= 𝑔𝐞𝖳{3},3 𝐑(𝑡𝑘 )𝐞{3},3 ∈ ℝ. (6) 5. Observability analysis

Because the accelerometer bias is assumed constant, we take In this section, we will study the observability of the
in consideration the invariance DT-LTV system (13) in a finite-time interval denoted by
𝑘0 ,𝓁 = [𝑘0 , 𝑘0 + 1, … , 𝑘0 + 𝓁], for 𝓁 ∈ ℕ.
𝑏[𝑘 + 1] = 𝑏[𝑘]. (7) Our study is divided in two parts. In the first one, we shall
According to (5), let us consider now two sets of discrete- briefly demonstrate the existence of specific necessary con-
time time-invariant auxiliary states given by ditions to guarantee observability when (17) in Assumption
( ) 1 is not satisfied. In the second, more general part, we will
𝑥1,𝑖 [𝑘] ∶= 𝛽𝑖 cos 𝜙𝑖 ∈ ℝ, and (8) show that (17) in Assumption 1 is a sufficient condition to
( ) attain observability.
𝑥2,𝑖 [𝑘] ∶= 𝛽𝑖 sin 𝜙𝑖 ∈ ℝ, (9)
The DT-LTV system (13) is of order 2𝑁 + 1, and, for
such that, based on (3) and (4), we may write, for 𝑖 = 𝓁 ≥ 2𝑁, the observability matrix associated with the pair
1, 2, … , 𝑁 (𝐀, 𝐂[𝑘]) can be written as (cf. Rugh (1996))

∑
𝑁
( ) ⎡ 𝐂[𝑘0 ] ⎤
𝑎𝑚 [𝑘] − 𝑢[𝑘] = sin 𝜔𝑖 𝑘𝑇𝑠 𝑥1,𝑖 [𝑘] ⎢ 𝐂[𝑘 + 1]𝐀 ⎥
𝑖=1 ⎢ 0 ⎥
(10) O ∶= ⎢ 𝐂[𝑘0 + 2]𝐀2 ⎥ ∈ ℝ(𝓁+1)×(2𝑁+1) . (18)
∑
𝑁
( ) ⎢ ⋮ ⎥
+ cos 𝜔𝑖 𝑘𝑇𝑠 𝑥2,𝑖 [𝑘]. ⎢𝐂[𝑘 + 𝓁]𝐀𝓁 ⎥
𝑖=1 ⎣ 0 ⎦

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 3 of 9

 Discrete-time Kalman lter for heave motion estimation
If the DT-LTV system (13) is observable, then O in (18) Geometrically, (22) establishes quite an interesting result.
must be full rank, in the sense that the initial condition 𝐱[𝑘0 ] It states that the sum of a sequence of real inner products
is uniquely determined from the collected set of measure- between rotating vectors and constant vectors related to 𝐝1
ments given by {𝐲[𝑘0 ], 𝐲[𝑘0 + 1], … , 𝐲[𝑘0 + 𝓁]}. and 𝐝2 must be always equal to −𝑑3 . Alternatively, it can
We begin by analyzing the case of a single sinusoidal be seen as a time-invariant inner product in a real vector
(𝑁 = 1, and thus, in view of (11), 𝝎 ≡ 𝜔 = 2𝜋𝑓 ), for which space. In order to better infer this last result, let us define
the minimum interval of observability, 𝑘0 ,2 , contains three an auxiliary constant vector given by
samples, meaning that, according to (18), we have [ ( )]
sin (𝝎𝑘0 𝑇𝑠 )
( ) ( ) 𝐯= ∈ ℝ2𝑁 . (23)
⎡ sin cos 𝝎𝑘0 𝑇𝑠
( 𝜔𝑘0 𝑇𝑠 ) ( 𝜔𝑘0 𝑇𝑠 ) 1⎤
cos
O = ⎢sin(𝜔(𝑘0 + 1)𝑇𝑠 ) cos(𝜔(𝑘0 + 1)𝑇𝑠 ) 1⎥ . (19) For each sample (whose instance of measurement is encoded
⎢ ⎥
⎣sin 𝜔(𝑘0 + 2)𝑇𝑠 cos 𝜔(𝑘0 + 2)𝑇𝑠 1⎦ by 𝓁), this constant vector 𝐯 is subject to a transformation
𝐐𝓁 ∈ ℝ2𝑁×2𝑁 given by
For O in (19) to be full rank, it suffices to show that [ ( ) ( )]
det(O) ≠ 0, for all 𝑘0 . Straightforward diag(cos( 𝝎𝓁𝑇𝑠) ) − diag(sin( 𝝎𝓁𝑇𝑠) )
( ) algebraic
( )manipu- 𝐐𝓁 ∶= . (24)
lations lead to det(O) = sin 2𝜔𝑇𝑠 − 2 sin 𝜔𝑇𝑠 , which diag(sin 𝝎𝓁𝑇𝑠 ) diag(cos 𝝎𝓁𝑇𝑠 )
is non-zero if and only if 𝜔𝑇𝑠 ≠ 𝑚𝜋, for 𝑚 ∈ ℕ∖{0} or, Therefore, using compact notation, we can rewrite (22) as
equivalently, if and only if 𝑓𝑠 ≠ 2𝑓 ∕𝑚, for 𝑚 ∈ ℕ∖{0}. This
outcome means that the sampling frequency can be smaller 𝐯𝖳 𝐐𝓁 𝐰 = −𝑑3 , for 𝓁 = 0, 1, 2, … , 2𝑁, (25)
than the actual wave frequency 𝑓 , as long as 𝑓𝑠 is not a where, for notational simplicity, we have defined
submultiple of 2𝑓 . In other words, (17) in Assumption 1 is [ ]
not a necessary condition to attain observability, but if it is 𝐝
𝐰 ∶= 1 ∈ ℝ2𝑁 . (26)
satisfied, in which case we have 𝑓𝑠 > 2𝑓 , then the DT-LTV 𝐝2
system (13) is observable for 𝑁 = 1 given any 𝑘0 ,2 . Note that 𝐐𝓁 is an orthonormal matrix in a real vector space,
For 𝑁 = 2, with 𝝎 ≡ [𝜔1 𝜔2 ]𝖳 = 2𝜋[𝑓1 𝑓2 ]𝖳 , meaning that ‖𝐐𝓁 ‖ = 1 for all 𝓁. Moreover, for 𝓁 = 0, we
computing det(O) by hand is still tractable; tedious but have 𝐐0 = 𝐈2𝑁 . For consequent values of 𝓁, and in order
straightforward algebraic manipulations lead to to guarantee (25) is satisfied, 𝐐𝓁 must either be an identity
( ) ( )[ ( ) ] matrix or, in the case when 𝑑3 = 0, meaning that 𝐯 and 𝐰
det(O) = −16 sin 𝜔1 𝑇𝑠 sin 𝜔2 𝑇𝑠 cos 𝜔1 𝑇𝑠 − 1 are orthogonal, be a reflection matrix. These two situations
[ ( ) ( )]2 [ ( ) ] (20)
cos 𝜔1 𝑇𝑠 − cos 𝜔2 𝑇𝑠 cos 𝜔2 𝑇𝑠 − 1 . correspond to rotations by angles of 2𝜋 and 𝜋, respectively.
However, if (17) in Assumption 1 is satisfied, they never
Based on (20), we conclude that if det(O) ≠ 0, then three occur. Hence, we conclude that (17) is a sufficient condition
conditions must be satisfied simultaneously, for 𝑚 ∈ ℕ∖{0}: to ensure the DT-LTV system (13) is observable for all 𝑘0 ,𝓁 .
i) 𝑓𝑠 ≠ 𝑓1 ∕𝑚; ii) 𝑓𝑠 ≠ 𝑓2 ∕𝑚; and iii) 𝑓𝑠 ≠ (𝑓2 −𝑓1 )∕𝑚. Once
again, (17) in Assumption 1 is not a necessary condition
to attain observability, but if it is satisfied, meaning 𝑓𝑠 > 6. DT-LTV Kalman filter design
2𝑓2 > 2𝑓1 , then the DT-LTV system (13) is observable In this section, we present the details of a Kalman filter
when 𝑁 = 2, for any 𝑘0 ,4 . implementation as a natural estimation solution for the DT-
Unfortunately, for 𝑁 ≥ 3, obtaining a closed form LTV system (13). We first recall the nominal measurements
analytic expression for det(O) becomes impracticable. A equation established in (4). In practice, the accelerometer
more clear-cut alternative to demonstrate observability lies readings are corrupted by noise, which we shall denote here
in a proof by contraposition, which we proceed to elaborate. by 𝑛𝑎 (𝑡𝑘 ) ∈ ℝ. This variable is assumed to represent a zero-
If O is not full rank, then there must exist a unit vector mean white Gaussian sensor noise sequence.
𝐝 ∈ ℝ2𝑁+1 , ‖𝐝‖ = 1, with 𝐝 ∶= [𝐝1 , 𝐝2 , 𝑑3 ], where Remark 1. The term 𝐞𝖳{3},3 𝐑(𝑡𝑘 )𝐞{3},3 in (4), which en-
𝐝1 , 𝐝2 ∈ ℝ𝑁 and 𝑑3 ∈ ℝ, such that O𝐝 = 𝟎(𝓁+1)×1 .
codes the pitch and roll angles of the platform, represents
Hence, according to (18), if the DT-LTV system (13) is not
an additional sensor measurement and, therefore, is also
observable, then, by using (14), it must be
corrupted by noise in practice. For filter design purposes, we
⎡ 𝐂[𝑘0 ]𝐝 ⎤ ⎡ 𝐂[𝑘0 ]𝐝 ⎤ ⎡0⎤ assume that the noise associated with this term is contained
⎢ 𝐂[𝑘 + 1]𝐀𝐝 ⎥ ⎢ 𝐂[𝑘 + 1]𝐝 ⎥ ⎢0⎥ in 𝑛𝑎 (𝑡𝑘 ). Admittedly, this design choice misrepresents the
⎢ 0
2 ⎥
0
underlying assumption by which the Kalman filter errors
⎢ 𝐂[𝑘0 + 2]𝐀 𝐝 ⎥ = ⎢⎢ 𝐂[𝑘0 + 2]𝐝 ⎥⎥ = ⎢⎢0⎥⎥ . (21)
follow a Gaussian distribution and, consequently, overrules
⎢ ⋮ ⎥ ⋮ ⋮
⎢𝐂[𝑘 + 𝓁]𝐀𝓁 𝐝⎥ ⎢⎣𝐂[𝑘 + 𝓁]𝐝⎥⎦ ⎢⎣0⎥⎦ claims of optimality.
⎣ 0 ⎦ 0
Similar to (12), let the comprehensive system state esti-
From the previous result, and in view of (16), we first obtain, mate be denoted by ̂ 𝐱[𝑘] ∈ ℝ2𝑁+1 , with
for 𝓁 = 0, 1, 2, … , 2𝑁, and for all 𝑘0 , [
̂
𝐱[𝑘] ∶= 𝑥 ̂1,1 [𝑘] 𝑥
̂1,2 [𝑘] … 𝑥
̂1,𝑁 [𝑘]
( )𝖳 ( )𝖳 ]𝖳 (27)
sin 𝝎(𝑘0 + 𝓁)𝑇𝑠 𝐝1 +cos 𝝎(𝑘0 + 𝓁)𝑇𝑠 𝐝2 = −𝑑3 . (22) ̂2,1 [𝑘] 𝑥
𝑥 ̂2,𝑁 [𝑘] ̂
̂2,2 [𝑘] … 𝑥 𝑏[𝑘] .

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 4 of 9

 Discrete-time Kalman lter for heave motion estimation
The discrete-time equations of a classic Kalman filter for Table 2
the DT-LTV system (13) are thus given by Wave frequency, amplitude and phase
( ) 𝑖
𝜔𝑖
(Hz) 𝐴𝑖 (m) 𝜙𝑖 (rad)
⎧̂𝐱[𝑘+1] = 𝐀̂𝐱[𝑘] + K[𝑘] 𝐲[𝑘] − 𝐂[𝑘]̂
𝐱[𝑘] , (28a) 2𝜋

⎪ ( 𝖳
) 𝖳
1 0.05 0.6853 -0.1834
⎪ K[𝑘] = 𝐀𝐏[𝑘]𝐀 + 𝐐 𝐂[𝑘] × 2 0.1 0.9182 -0.0421
⎨ ( ( ) )−1 3 0.15 0.1323 1.1805
⎪ × R + 𝐂[𝑘] 𝐀𝐏[𝑘]𝐀𝖳 + 𝐐 𝐂[𝑘]𝖳 , (28b)
⎪ ( )( ) 4 0.2 0.5896 -2.9270
𝖳
⎩𝐏[𝑘+1] = 𝐈2𝑁+1 − K[𝑘]𝐂[𝑘] 𝐀𝐏[𝑘]𝐀 + 𝐐 , (28c) 5 0.25 1.3785 -2.5686
where Q ∈ ℝ(2𝑁+1)×(2𝑁+1) , Q ≻ 𝟎, and R > 0, are the co-
variances of the process and observation noises, respectively. where 𝑛𝑢𝑝 (𝑡), 𝑛𝑢𝑟 (𝑡) ∈ ℝ represent another two sequences of
Each of these two covariances, herein assumed constant, is additive white Gaussian noise with standard deviation of 1
associated with a different additive white Gaussian noise degree. The evolution of 𝑢(𝑡𝑘 ) is showcased in Fig. 1.
distribution, and can be seen as a tuning knob adjusted for
10
the best performance.
Finally, an estimate of the heave motion follows as 8

Gravitational acceleration m/s2

6
∑𝑁
1 { ( ) ( ) }
𝑧̂[𝑘] = − sin 𝜔𝑖 𝑘𝑇𝑠 𝑥̂1,𝑖 [𝑘] + cos 𝜔𝑖 𝑘𝑇𝑠 𝑥̂2,𝑖 [𝑘] . 4
2
𝑖=1 𝜔𝑖
2
(29)
0

7. Simulation Validation and Testing -2

7.1. Case 1 - Discrete wave spectrum with known -4

0 5 10 15 20 25 30
frequencies Time (s)

We set out to validate our strategy, particularly by as- Figure 1: Evolution of 𝑢(𝑡𝑘 ) according to (30) and (31).
sessing the performance of the Kalman filter (28) in the
presence of sensor noise, including cross-correlation terms,
as stressed in Remark 1. In regard to the Kalman filter settings, the initial state
Consider a mission scenario where a surface vessel estimate is set to zero, i.e., ̂
𝐱(𝑡𝑘0 = 0) ∶= 𝟎2𝑁+1 . The
equipped with, e.g., a winch system, is subjected to an initial error covariance matrix, and the constant process and
oscillatory motion consisting in the overlay of five waves observations covariances are set (for 𝑁 = 5) as
whose properties are listed in Table 2. Further suppose
𝐏(𝑡𝑘0 = 0) ∶= blkdiag(100𝐈𝑁 , 100𝐈𝑁 , 10−1 ),
that the vessel is equipped with an IMU whose readings
of acceleration are corrupted by additive white Gaussian 𝐐 ∶= blkdiag(10−10 𝐈𝑁 , 10−10 𝐈𝑁 , 10−5 ), (32)
noise with standard deviation of 0.01 m/s2 , as well as by a 𝐑 ∶= 0.01,
sensor offset characterized√by a Brownian motion governed
by 𝑏(𝑡𝑘+1 ) = 𝑏(𝑡𝑘 ) + 0.003 𝑇𝑠 𝑟𝑏 (𝑡𝑘 ), with 𝑏(𝑡𝑘0 = 0) = 0.05 respectively. Finally, the sampling rate is set as 𝑓𝑠 = 50
m/s2 and 𝑟𝑏 ∈ ℝ being a sequence of normally distributed Hz. Note that the matrices in (32) were empirically adjusted
random numbers. We assume that the motion of the sea for best performance, but they still reflect the knowledge of
waves causes the surface vessel to pitch and roll as described model and sensor specifications. In particular, the third entry
by the following parametric functions of time: of 𝐐 is set larger than the remaining ones to allow the filter
( ) to track the time-varying nature of 𝑏(𝑡𝑘 ).
𝜃(𝑡𝑘 ) = 0.92 sin 2𝜋0.033𝑡𝑘 + 2𝜋(𝑟𝑝1 (𝑡𝑘 ) − 0.5) Let us first present the results associated with the esti-
( ) mation error regarding 𝑥1,𝑖 and 𝑥2,𝑖 , for 𝑖 = 1, 2, … , 𝑁,
+ 0.15 sin 2𝜋0.165𝑡𝑘 + 2𝜋(𝑟𝑝2 (𝑡𝑘 ) − 0.5)
( ) plotted in Fig. 2. The evolution of the estimation error
+ 0.53 sin 2𝜋0.32𝑡𝑘 + 2𝜋(𝑟𝑝3 (𝑡𝑘 ) − 0.5) and, associated with the heave is plotted in Fig. 3. Notice the
( ) (30)
𝜑(𝑡𝑘 ) = 0.18 cos 2𝜋0.066𝑡𝑘 + 2𝜋(𝑟𝑟1 (𝑡𝑘 ) − 0.5) filter’s quick convergence, with steady-state behavior being
( ) reached roughly around 10 seconds. The filter’s accuracy in
+ 1.15 cos 2𝜋0.18𝑡𝑘 + 2𝜋(𝑟𝑟2 (𝑡𝑘 ) − 0.5)
( ) steady-state is approximately 4 cm, which is a very good
+ 0.72 cos 2𝜋0.25𝑡𝑘 + 2𝜋(𝑟𝑟3 (𝑡𝑘 ) − 0.5) , result for this kind of application.
where 𝑟𝑝𝑖 , 𝑟𝑟𝑖 ∈ [0, 1], for 𝑖 = 1, 2, 3, are uniformly The estimation over time of each of the five wave am-
distributed random numbers, with 𝜃(𝑡) and 𝜑(𝑡) representing plitudes and five wave phases is showcased in Figs. 4 and 5,
the pitch and roll angles, respectively, such that, according respectively. In line with the results shown in Figs. 2 and 3,
to (6), we have once more the quick convergence is evidenced by the short
( ) ( ) transient period, which lasts roughly 10 seconds. Most in-
𝑢(𝑡𝑘 ) = 𝑔 cos 𝜃(𝑡) + 𝑛𝑢𝑝 (𝑡) cos 𝜑(𝑡) + 𝑛𝑢𝑟 (𝑡) , (31) terestingly, we observe that the frequency associated to each

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 5 of 9

 Discrete-time Kalman lter for heave motion estimation
40 100

20 50
0
0
-20
-50

Phase (deg)
-40
-100
40
-150
20
-200
0
-250
-20

-40 -300
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Time (s) Time (s)

Figure 2: Evolution of estimation errors 𝑥̃1,𝑖 (𝑡𝑘 ) ∶= 𝑥1,𝑖 (𝑡𝑘 ) − Figure 5: Evolution of the phase estimates.
̂1,𝑖 (𝑡𝑘 ) and 𝑥
𝑥 ̂2,𝑖 (𝑡𝑘 ), for 𝑖 = 1, 2, … , 𝑁 .
̃2,𝑖 (𝑡𝑘 ) ∶= 𝑥2,𝑖 (𝑡𝑘 ) − 𝑥
Dashed lines correspond to the evolution of the +∕− standard
deviation of the errors.
Note that the filter is able to track relatively fast changes with
15 an accuracy typically less than 10%.
40
10

30
5
20
Error (m)

Amplitude (mg)
0 10

0
-5 0.04
0.02
0
-0.02 -10
-0.04 7.5
-10 20 22 24 26 28 30
-20 7

-15 6.5
-30
0 5 10 15 20 25 30
Time (s) 6
20 22 24 26 28 30
-40
0 5 10 15 20 25 30
Figure 3: Evolution of heave error 𝑧̃(𝑡𝑘 ) ∶= 𝑧(𝑡𝑘 ) − 𝑧̂(𝑡𝑘 ). Time (s)

Figure 6: Evolution of bias estimate ̂𝑏(𝑡𝑘 ), alongside true 𝑏.

amplitude, and not amplitude’s actual value, impacts on the

transient performance: each wave amplitude is characterized
Remark 2. The interval between frequencies has an impact
by larger deviations than the succeeding one.
on the performance of the filter. More specifically, narrower
18 gaps lead to slower convergence for the same sampling rate.
16 By increasing the sampling rate it is possible to mitigate
14
that effect. However, the impact caused by narrower gaps
12
1.5 between frequencies is not too critical in terms of relative
change in standard deviation of the heave estimation error.
Amplitude

10 1

8 To corroborate the claim in the previous remark, let

0.5
6 us run our filter 100 times, where, for each simulation,
4 0 we consider only two wave frequencies that are distanced
20 22 24 26 28 30
2 from 0.0005 Hz up to 0.05 Hz in steps of 0.0005 Hz. The
0
parameters across all simulations are identical to the ones
0 5 10 15 20 25 30
Time (s)
reported in this section, including the time-varying bias drift
governed by Brownian motion. In Fig. 7, we present the
Figure 4: Evolution of the amplitude estimates. mean and standard deviation, computed in steady state, for
15 < 𝑡 < 30 seconds, in function of Δ𝑓 . The trend is evident:
Finally, the bias estimate is shown in Fig. 6, alongside for the same sampling rate, higher gap between frequencies
the actual sensor offset considered in the simulation. For lead to higher accuracy within the same period of time, but,
convenience of representation, the units are displayed in mg. overall, the filter performance is satisfactory.
The dashed lines mark the +∕− standard deviation from
the current estimate, as obtained from the evolution of the
7.2. Case 2 - Continuous wave spectrum
To further demonstrate the usefulness of the proposed
Kalman filter’s covariance matrix 𝐏’s last diagonal entry.
methodology, especially the fact that it does not require any
additional computations related to the Fast fourier transform,

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 6 of 9

 Discrete-time Kalman lter for heave motion estimation
15
Table 3
10
Case 2: Heave and bias root mean square error
N Heave (m) Bias (mg)
5 50 0.8525 1.7161
55 0.7344 1.4620
0 60 0.8423 1.6734
65 0.9611 1.7490
-5 70 0.7772 1.3876
Mean (cm) 75 0.5393 1.0524
Standard deviation (cm)
-10
80 0.6767 1.1725
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05

Figure 7: Steady-state mean and standard deviation vs. gap Assumption 2. The 𝑁 wave frequencies comprise most
between frequencies 𝑓1 = 0.05 Hz and 𝑓2 = 𝑓1 + Δ𝑓 Hz. of the energy of the continuous wave spectrum. Therefore,
according to (1) and (5), we establish the approximation:

consider now a mission scenario where the oscillatory mo- ∑𝑁

tion of the sea waves, characterized by a Sea State Code 1 { ( ) ( ) }
𝑧[𝑘] ≈ − sin 𝜔𝑖 𝑘𝑇𝑠 𝑥1,𝑖 [𝑘] + cos 𝜔𝑖 𝑘𝑇𝑠 𝑥2,𝑖 [𝑘] .
2
7, is described by a continuous Pierson-Moskowitz wave 𝑖=1 𝜔𝑖
spectrum given by (cf. Fossen (2011)[Chapter 8]) (34)
( )
8.1 × 10−3 𝑔 2 3.11 1 Accordingly, suppose 𝜔1 and 𝜔𝑁 are the smallest and
𝑆(𝑓 ) = exp − , (33) largest frequencies of the interval of interest, respectively.
(2𝜋𝑓 )5 𝐻𝑠2 (2𝜋𝑓 )4
The corresponding frequency step size is thus given by
where 𝐻𝑠 ∈ [6, 9] m is the significant wave height for Sea 𝜔𝑁 − 𝜔1
State Code 7. A continuous spectrum of this type typically Δ𝜔 ∶= , (35)
𝑁 −1
resembles the one shown in Fig. 8.
which allows us to slightly redefine the vector of frequencies,
1
first introduced in (11), as a vector of equally spaced values
0.9
given by
0.8
[ ]
Normalized Amplitude

0.7
𝝎 ∶= 𝜔1 𝜔1 + Δ𝜔 𝜔1 + 2Δ𝜔 … 𝜔𝑁 ∈ ℝ𝑁 . (36)
0.6

0.5
Note that this uniform spacing of frequencies reflects the
0.4
fact that no specific knowledge about the wave spectrum is
0.3
obtained a priori.
0.2
We then update the covariance of the process noise in the
0.1
Kalman filter as 𝐐 ∶= blkdiag(10−2 𝐈𝑁 , 10−2 𝐈𝑁 , 10−2 ), and
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 we set 𝜔1 = 0.1𝜋 and 𝜔𝑁 = 0.8𝜋. Sensor specifications and
f (Hz)
other simulation parameters remain unchanged.
Figure 8: A simulated example of a Pierson-Moskowitz spec- In Fig. 9 we present the evolution of the resulting heave
trum of ocean waves described by Sea State Code 7. motion associated with (33) and the corresponding ship
dynamics, as well as the evolution of several estimation
Consider as well a large vessel, namely the HMS Nor- sequences obtained for different values of 𝑁 > 50. The
folk, weighing 2500 tonnes, of dimensions 137 × 15 × 16 evolution of the sensor bias and corresponding estimation
m, subjected to the aforementioned waves. The incidence sequences is displayed in Fig. 10.
angle, i.e., the angle between the waves and the direction of As expected, a larger 𝑁 yields more accurate results
motion of the ship was set to 60 degrees. The dynamics of since the whole continuous spectrum is better approximated.
this surface vessel was computed through the Ship Simulator Nevertheless, the overall Kalman filter performance is sat-
available in Granström and Bernat (2021) to obtain realistic isfactory, as evidenced by the calculated root mean square
sequences of pitch and roll, as well as heave motion data. errors listed in the Table 3. Obtaining a similar performance
The region of the spectrum that is of interest is assumed for 𝑁 < 50 is also possible to some extent, although
relatively well-known a priori. In view of the discrete- technically more difficult since the approximation in (34)
time nature of our strategy, we consider that region as an grows weaker. Further tweaking of the initial covariance
interval of 𝑁 discrete frequencies, where, as opposed to the matrix 𝐏, as well as of the constant covariance 𝐑 in (28)
scenario described in Section 7.1, the value of 𝑁 can be may help but with limited efficacy.
arbitrarily chosen relying on an educated guess. This design
trait motivates the following assumption.

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 7 of 9

 Discrete-time Kalman lter for heave motion estimation
10
We deliberately selected an insufficient range to high-
light these shortcomings, which are reflected in the smaller
5
accuracy listed in Table 3, which compares slightly poorer
than the strictly discrete spectrum case, as seen from Fig. 3.
Amplitude (m)

-5 0.9

0.8

Normalized Amplitude
-10 0.7

0.6

-15 0.5
0 50 100 150 200 250 300 350
Time (s) 0.4

Figure 9: Simulated heave based on JONSWAP spectrum and 0.3

corresponding estimation sequences for dierent values of 𝑁 . 0.2

0.1

25 0
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
f (Hz)
20
Figure 11: Continuous Pierson-Moskowitz spectrum alongside
reconstructed discrete spectra based on estimated amplitudes
obtained for dierent values of 𝑁 .
15
Amplitude (mg)

10

8. Conclusions
0
This paper presented a discrete-time linear time-varying
-5
0 50 100 150 200 250 300 350
Kalman filter to solve the problem of heave motion esti-
Time (s) mation using biased accelerometer measurements. Making
Figure 10: Simulated bias drift with Brownian motion and use of convenient and well-known properties of the sin(∙)
corresponding estimation sequences for dierent values of 𝑁 . function, we developed a general framework encapsulating
a linear system whose state encompasses implicit repre-
sentations of wave amplitudes and phase shifts, as well
For the sake of completeness, we display in Fig. 11 a as a constant sensor bias. We showed that the system is
discrete-time normalized reconstruction of the signal spec- observable, which guarantees the Kalman-based estimator
trum for each value of 𝑁. As expected, these discrete se- is globally asymptotically stable. The proposed strategy was
quences closely resemble the continuous spectrum, origi- validated with realistic numerical examples, including an
nally plotted in Fig. 8, apart from some excited modes related accurate representation of a continuous Pierson-Moskowitz
to the motion of the ship. wave spectrum within an ocean context where a large surface
We emphasize the fact that the choice of frequency vessel is subjected to wave motion. Longer missions, as
interval, i.e., 𝜔1 and 𝜔𝑁 , can have an impact on the filter’s well as greater sailed distances, will certainly see the sea
estimation performance, most evidently inferred through the state changing over time, which will cause different spectra
gap in Fig. 11, particularly the empty area to the left of to manifest. Adapting our proposed methodology to this
𝑓 = 0.05 Hz. Allowing 𝝎 to encompass frequencies too scenario will be considered in future work.
close to zero will lead to some ambiguity in the estimation,
as the filter struggles to separate data from low frequencies References
and data from the sensor drift, which is caused by an ill- Auestad, Ø.F., Gravdahl, J.T., Fossen, T.I., 2013. Heave Motion Estimation
conditioned observability matrix. on a Craft Using a Strapdown Inertial Measurement Unit. IFAC Proceed-
Moreover, we stress that the proposed strategy consists ings Volumes 46, 298–303. doi:10.3182/20130918-4-jp-3022.00033.
in a fully discrete-time implementation and that, according Batista, P., Oliveira, P., Silvestre, C., 2018. A Globally Exponentially
to (34), a finite number of frequencies can approximate a Stable Solution for Frequency Estimation, in: 2018 IEEE Conference
on Decision and Control (CDC), IEEE. doi:10.1109/cdc.2018.8619093.
continuous spectrum, but it only does so well when said Ben, Y., Gao, Q., Wei, T., Gong, S., Li, Q., 2022. Real-time heave motion
number is relatively large, and the range defined by 𝜔1 measurement by adaptive band-pass filter based on strapdown INS.
and 𝜔𝑁 encloses the whole spectrum. The frequencies of Ocean Engineering 262, 112278. doi:10.1016/j.oceaneng.2022.112278.
principal excited modes do not, in most cases, coincide with Blake, S.J., 2007. Heave compensation using time-differenced carrier
the selected frequencies, in which case the Kalman filter observations from low cost GPS receivers. Ph.D. thesis. University of
Nottingham.
benefits from of its well adjusted process covariance in order Collins, C.O., Lund, B., Waseda, T., Graber, H.C., 2014. On recording sea
to track slowly-time varying changes: a wrong frequency is surface elevation with accelerometer buoys: lessons from ITOP (2010).
equivalent to an unknown time-varying phase. Ocean Dynamics 64, 895–904. doi:10.1007/s10236-014-0732-7.

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 8 of 9

 Discrete-time Kalman lter for heave motion estimation
El-Hawary, F., 1982. Compensation for source heave by use of a Kalman
filter. IEEE Journal of Oceanic Engineering 7, 89–96. doi:10.1109/joe.
1982.1145516.
Fossen, T.I., 2011. Handbook of Marine Craft Hydrodynamics and Motion
Control. Wiley. doi:10.1002/9781119994138.
Granström, L.H., Bernat, M., 2021. Ship Simulator. https://github.com/
LINK-SIC-2021-Bernat-Granstrom/ship-simulator.
Huang, B., Jiang, J., Zou, Z., 2021. Online Prediction of Ship Coupled
Heave-Pitch Motions in Irregular Waves Based on a Coarse-and-Fine
Tuning Fixed-Grid Wavelet Network. Journal of Marine Science and
Engineering 9, 989. doi:10.3390/jmse9090989.
Johansen, T.A., Fossen, T.I., Sagatun, S.I., Nielsen, F.G., 2003. Wave
synchronizing crane control during water entry in offshore moonpool
operations-experimental results. IEEE Journal of Oceanic Engineering
28, 720–728. doi:10.1109/joe.2003.819155.
Kielland, P., Hagglund, J., 2015. Using DGPS to Measure the Heave Motion
of Hydrographic Survey Vessels. The International Hydrographic
Review 72. URL: https://journals.lib.unb.ca/index.php/ihr/article/
view/23195.
Küchler, S., Eberharter, J.K., Langer, K., Schneider, K., Sawodny, O.,
2011a. Heave Motion Estimation of a Vessel Using Acceleration
Measurements. IFAC Proceedings Volumes 44, 14742–14747. doi:10.
3182/20110828-6-it-1002.01935.
Küchler, S., Mahl, T., Neupert, J., Schneider, K., Sawodny, O., 2011b.
Active Control for an Offshore Crane Using Prediction of the Vessel’s
Motion. IEEE/ASME Transactions on Mechatronics 16, 297–309.
doi:10.1109/tmech.2010.2041933.
Li, J., Li, L., Wen, X., 2022. A novel method of estimating ship heave
motion based on dual inertial measurement units. Journal of Physics:
Conference Series 2183, 012002. doi:10.1088/1742-6596/2183/1/012002.
Li, S., Wei, J., Guo, K., Zhu, W., 2017. Nonlinear Robust Prediction
Control of Hybrid Active–Passive Heave Compensator With Extended
Disturbance Observer. IEEE Transactions on Industrial Electronics 64,
6684–6694. doi:10.1109/tie.2017.2698358.
MacGregor, 2022. MacGregor Deck Handling Solutions. https://www.
macgregor.com/. Accessed: 2022-10-20.
Nielsen, U.D., Brodtkorb, A.H., Sørensen, A.J., 2019. Sea state estimation
using multiple ships simultaneously as sailing wave buoys. Applied
Ocean Research 83, 65–76. doi:10.1016/j.apor.2018.12.004.
Richter, M., Schaut, S., Walser, D., Schneider, K., Sawodny, O., 2017.
Experimental validation of an active heave compensation system: Esti-
mation, prediction and control. Control Engineering Practice 66, 1–12.
doi:10.1016/j.conengprac.2017.06.005.
Rugh, W.J., 1996. Linear system theory. Prentice Hall.
Skaare, B., Egeland, O., 2006. Parallel Force/Position Crane Control in
Marine Operations. IEEE Journal of Oceanic Engineering 31, 599–613.
doi:10.1109/joe.2006.880394.
Tang, G., Yao, X., Li, F., Wang, Y., Hu, X., 2022. Prediction about
the vessel’s heave motion under different sea states based on hybrid
PSO_ARMA model. Ocean Engineering 263, 112247. doi:10.1016/
j.oceaneng.2022.112247.
Wago, 2022. Wago Deck Handling Cranes. https://www.wago.com/global/
deck-handling-cranes. Accessed: 2022-10-20.
Woodacre, J.K., Bauer, R.J., Irani, R.A., 2015. A review of vertical
motion heave compensation systems. Ocean Engineering 104, 140–154.
doi:10.1016/j.oceaneng.2015.05.004.
Xie, Y., Zhao, X., Luo, M., 2022. An active-controlled heaving plate break-
water trained by an intelligent framework based on deep reinforcement
learning. Ocean Engineering 244, 110357. doi:10.1016/j.oceaneng.
2021.110357.
Zhou, H., Cao, J., Yao, B., Lian, L., 2021. Hierarchical NMPC–ISMC
of active heave motion compensation system for TMS–ROV recovery.
Ocean Engineering 239, 109834. doi:10.1016/j.oceaneng.2021.109834.

Reis, Batista, Oliveira and Silvestre: Preprint submitted to Elsevier Page 9 of 9