# Reconstruction of Riser Profiles by an Underwater Robot Using Inertial Navigation

Authors: Luciano Luporini Menegaldo; Stefano Panzieri; Cassiano Neves
Journal: Journal of Robotics
 (2012)

---

## Abstract

This paper proposes a kinematic model and an inertial localization system architecture
for a riser inspecting robot. The robot scrolls outside the catenary riser, used for underwater
petroleum exploration, and is designed to perform several nondestructive tests. It
can also be used to reconstruct the riser profile. Here, a realistic simulation model of robot
kinematics and its environment is proposed, using different sources of data: oil platform
characteristics, riser static configuration, sea currents and waves, vortex-induced vibrations,
and instrumentation model. A dynamic finite element model of the riser generates
a nominal riser profile. When the robot kinematic model virtually scrolls the simulated
riser profile, a robot kinematic pattern is calculated. This pattern feeds error models of a
strapdown inertial measurement unit (IMU) and of a depth sensor. A Kalman filter fuses
the simulated accelerometers data with simulated external measurements. Along the riser
vertical part, the estimated localization error between the simulated nominal and Kalman
filter reconstructed robot paths was about 2 m. When the robot model approaches the
seabed it assumes a more horizontal trajectory and the localization error increases significantly.


---

## Body

## 1.  Introduction

One of the key elements of deep-water petroleum exploration is the production riser. Risers are the ducts that transport petroleum, water or gases from the exploitation well up to the production platform. Either rigid or flexible types of risers may be used in the oil field. Both types are submitted to a broad spectrum of failure causes [1]: mechanical loads, aging, corrosion, erosion, temperature effects, installation or fabrication nonconformities, and so forth. Therefore, the availability of inspection tools to assess riser integrity status in situ is highly desirable. Such procedures are performed mainly by visual inspection with remotely operated vehicles (ROVs) [2] or autonomous underwater vehicles (AUVs) [3]. In some cases, sensors are installed directly on fixed points of the riser surface to measure strain and riser motion [4]. Other types of nondestructive testing (NDT) techniques can be used, such as magnetic, radiographic, or ultrasound methods [5]. In these cases, however, the operational constraints for using human operators are a major problem. A few papers address robotic devices specifically designed for underwater riser inspection. Psarros and his collaborators [6] proposed a robot that moves along the riser by using a mechanism composed of two parts. One part stays attached to the riser body, and the other part moves towards the riser’s side, in a cyclical manner.A major technical problem in robotic underwater inspection is the navigation and/or localization of the robot in a highly dynamic sea environment. Navigation is especially critical for AUVs and somewhat critical for ROVs. Lee at al. [7] addressed this problem by using several sensors fused by a multirate Extended Kalman Filter (EKF). The sensors set included a strapdown inertial platform, a Doppler velocity log (DVL), magnetic compass, and a depth sensor. However, they had sonar transducers installed in an underwater reference station and in the remote vehicle. Jouffroy and Opderbecke [8] addressed the problem of measuring the horizontal position of a ROV by using a gyro-Doppler together with an ultrashort baseline (USBL) acoustic positioning system. Diffusion-based observers were used to process a trajectory segment, instead of a typical point-by-point localization. He et al. [9] proposed an approach based on an invariance extended Kalman filter (IEKF) to address the problems of using sonar in shallow waters. In the case when the robot is mechanically linked to the inspected structure, the key problem is to localize precisely where it is at every instant of time. Such localization coordinates are associated to NDTs data and the flaws position can be precisely determined.Recently, our group designed and built a prototype of a robotic device specifically designed to perform nondestructive testing (NDT) in production risers [10]. The robot has neutral floatability and embraces the riser by moving along its outside (Figure 1), using a pair of thrusters for propulsion as well as polymeric wheels to guarantee sliding and correct alignment with the riser surface. In Figure 2 it is possible to observe how the robot attaches the riser by opening and closing its motorized arms. This operation is assisted by a human diver. It communicates with the operator’s computer by means of an umbilical cable that transmits power, images and control commands. The dimensions of the robot and additional parameters used along the work are shown in Table 2. This robot will be able to perform several NDT procedures, such as ultrasound, imaging, and mechanical vibration measurements.Figure 1
Robot prototype.Figure 2
Transversal view of the robot frame showing how it attaches the riser, by opening and closing its motorized arms. The polymeric free wheels that effectively touch the riser surface are also shown.This paper proposes a kinematic model of the robot performing a riser profile cast mission, in a realistic simulated environment. Initially, a riser dynamic profile is estimated using a finite element model of the riser subjected to sea and ship motions. The nominal robot kinematic path (including position, velocity and acceleration), as it scrolls by the riser, iscontaminated with experimental errors, simulated by IMU and depth sensor models. The simulated sensor data is used by a Kalman filter to estimate the original robot path. This path is a good estimate of the actual riser profile, if robot mission time is small, compared to platform motion.The obtained profile can be used as an imposed displacement data for some structural analysis software based on finite elements techniques, that allows stress to be calculated. In addition, the localization algorithm can be used to associate each NDT measurement with its riser geometrical coordinate. These two aspects are intimately connected, and the localization algorithm can be used either to cast the profile, for fast robot runs, or localize the NDTs.Reproducing the expected environmental conditions, to test the proposed approach, in a laboratory experiment is essentially impracticable. Field tests, by his turn, should require a expensive positioning system such as a 3D sonar, which does not operate at the required frequency resolution, due to the presence of vortex induced vibrations (VIVs). Therefore, a simulation of the riser application, together with simulated sensors, was used to assess the performance of the localization algorithm.Actually, a particular environmental and riser configuration scenario is being addressed in this paper. However, the approach is likely to be applicable to similar situations. Other devices that move along subsea pipe systems, such as flowlines, jumpers, and umbilicals, could employ some of the main ideas presented in this paper. No additional localization devices, such as sonar beacons, are needed.The localization problem formulation used a standard Kalman filter as a sensor fusion algorithm based on a simple kinematic model of a strapdown IMU fused with a depth sensor. More sophisticated sensor fusion algorithms or state-space models for the system (e.g., dynamical models) could also be tested in future implementations, but the problem formulation would be probably quite similar.


## 2.  Riser Simulation Conditions

The particular sea and ship motion conditions selected for running the simulations corresponded to asevere condition, relatively similar to that typically found in Campos Basin, in the southeast of South America coast. Actually, they were designed to be worse than the most severe scenario in which a robot is expected to operate (Table 1). Additionally, under milder sea conditions, the localization performance should be expectable to be better than shown here.Table 1
Parameters used in the simulations.
ParameterValueRobot length1133 mmRobot max. outside diameter800 mmRobot mass73 kgRiser outer diameter295.5 mmRiser inner diameter203.5 mmRiser length1530 mRiser internal pressure60 barOuter riser drag coefficient (Cd)1.2Water depth1180 mSeabed axial friction0.35Seabed transversal friction0.9Seabed vertical stiffness104.3N/m2Wave height5 mWave period10 secondsFPSO length330 mTable 2
Identified parameters for the IMU sensor model1.
Variableτr(Rxx(s))τr(Allan(s))Low cut (Hz)High cut (Hz)σbiasσbwRoll°62920*5200.16350.0072Pitch°57820*6200.15530.0072Yaw°140320*5120.21440.0095AccX g0.644*100.335.2768e-40.0041AccY g0.581*100.2530.00110.0037AccZ g0.645*100.336.4654e-40.0033AngRateX°/s0.767*15——0.00760.0076AngRateY°/s0.754*20——0.00740.0074AngRateZ°/s0.745*10——0.00570.00571Rxx means that used time constant is from autocorrelation, and Allan from Allan variance plot. Low cut and high cut correspond to the cutting-off frequencies of the low-pass and high-pass filter used to calculate σbias and σbw, respectively. σbw and σbias have the same unities of their respective measurement variables.*Used in the simulations.Data from a flexible free-hanging riser installed in a PETROBRAS (Brazilian State oil company) Turret Floating Production Storage Offloading (FPSO) oil platform, which is currently in operation in the Campos Basin, was used as the inputs for the riser simulation software FLEXCOM. This is a finite elements software customized for nonlinear static and dynamic analysis of offshore systems, used worldwide by the petroleum industry from the last 20 years, and validated against experimental tests and other finite elements packages [11, 12]. The software allows riser responses to be simulated with several kinds of platform characteristics, sea current profiles, hydrodynamic loads, regular and irregular waves, and so forth. Two situations were studied: static and dynamic. In the static analysis, only the equilibrium configuration of the riser is considered, without motions other than from the robot itself. Three positions of the platform with relation to the wheel head were considered: standard and with a ship offset of 150 m in the directions near and far. For the dynamic analysis, sea conditions with regular waves, reaching the ship 45° obliquely, were used to generate the riser motion profiles.To estimate numerically the localization system performance and calculate the associated displacement errors, the arrangement shown in Figure3 was used. By simulating the kinematic model (Section 4), a physical profile of Euler angles, inertial accelerations, and depth was calculated. These simulated variables are then contaminated by noise, using the IMU and depth sensor instrumentation error models (Section 5), providing realistic sensor outputs. The trajectory run by the robot is estimated by a sensor fusion algorithm (Kalman filter), using such noisy sensor data (Section 6). Finally, both physical and estimated riser shapes are compared, to estimate the localization error along the run (Sections 7 and 8).Figure 3
Overview of the system simulation analysis. A dynamic riser profile is generated through an FEM analysis. Using robot’s kinematical and instrumentation model, the expected sensor readings are used to reconstruct the riser profile by the localization algorithm, which is compared to the original profile from FEM output.


## 3.  Architecture of the Localization System

The proposed localization system, that will be simulated numerically here, is shown in Figure4. An IMU measures three accelerations, angular velocities, and Euler angles from the robot as it scrolls along the riser. Accelerations are measured in the local reference frame (see definition in Section 4). Using the Euler angles, the accelerations are transformed to the global reference system, using the classic strapdown inertial navigation approach [13]. The simulated measurements from the IMU are fused by the Kalman Filter (KF) with the processed simulated external measurements from the depth sensor. Its output is an estimated state vector, that includes robot position, in the global reference frame.Figure 4
Architecture of the proposed localization system architecture. The robot scrolling along the riser generates a set of kinematic variables that is measured by the IMU. The accelerations are transformed to the global reference system, using IMU attitude outputs, and considered as the inputs of the KF. To compensate the drift caused by integrated sensor noise, the KF fuses the IMU with depth sensor data, that is, an absolute measurement. Before entering in the KF, the external measurements are processed by the method described in Section6, using IMU attitude data.


## 4.  Generation of Static and Dynamic Sensor Profiles

A simulated profile of inertial sensor physical excitation was obtained. First, thenominal or static riser geometry was obtained from the FEM (finite element model) analysis. The global reference system (XG, YG, ZG) was positioned in the turret center, such that X axis points to an arbitrary direction (e.g., north), Z points up, and Y is orthogonal to both (Figure 5). This reference frame is attached to the FPSO, being thus a slowly moving frame. However, due to the very low frequency of the FPSO movement compared to robot mission time span, the global reference frame was considered as being inertial.Side and upper views of the FPSO and riser configuration. Global (XG, YG, ZG) and local (Xl, Yl, Zl) reference systems position and orientation are shown. The FPSO and riser catenary plane are considered as being aligned with the sea current, forming an angle ϕcat with the XG direction in XGYG plane.
(a)(b)The local reference system (Xl, Yl, Zl) is localized in the geometric center of the riser section that coincides with the position, relative to robot height, where the inertial sensor is expected to be installed, and moves together with the robot. The center of the local reference system was defined, relatively to the global reference system, using the same nodes of the FEM mesh, which coordinates in the global reference system are PiG=(x(i),y(i),z(i)), i=1,…,N (N is number of nodes in the mesh and the superscript G express the reference system). Each Zl vector is found from a unit vector Du that connects node i to node i+1 of the FEM mesh. The negative sign for defining Zl in (2) is due to the fact that Du points down:(1)D=[x(i+1)-x(i)y(i+1)-y(i)z(i+1)-z(i)],Du=D‖D‖,(2)Zl=-Du.In the sequence, a particularYl was chosen (along with its orthogonal counterpart Xl), such that it could represent one plausible mission trajectory. In the present version of the model, the robot communicates with the operator in the platform through an umbilical cable, and only a small amount of spin can be allowed, to prevent the cable from curling along the riser. Therefore, Yl is simultaneously orthogonal to a point sufficiently far in the X direction and to −Zl. The Xl vector was orthogonal to both Zl and Yl. For each local reference system i, a directors cosine matrix (DCM) was defined as follows:(3)DCM=[X1xY1xZ1xX1yY1yZ1zX1zY1zZ1z],
where Xl=(X1x,X1y,X1z)T is the unit vector that defines Xl, and similarly Yl and Zl.Euler anglesψ, θ, and ϕ were defined between the XG and Xl, YG and Yl, ZG and Zl axis, respectively. In every step, these angles were derived from the DCM matrix, by using the formulas presented by [14]. The resulting transformation matrix from local to global reference system R1G=RxRyRz was posteriorly verified to be the same as DCM, in order to determine if singularities were present. A fixed point PS=(0,0.5,0) expressed in the local reference frame corresponds to the robot body point where the IMU will be possibly installed. This point was arbitrarily chosen within a reasonable Y distance from the riser section center, where the origin of the local reference frame is positioned. The coordinates of PS, expressed in the global coordinate frame, were assumed as the sensor displacement profile.


### 4.1. Sea Current Effect

A linear current profileVcur from the maximum sea current velocity (Vcurmax) at sea level to 0 at seabed [15] was considered as being aligned with the riser catenary plane. Vcurmax was assumed as 1.68 m/s, from Campos Basin data:(4)Vcur(i)=RZcat[0Vcurmax(z(N-1)-z(i))z(N-1)-z(1)0]i=1,…,N,
where RZcat is the rotation matrix associated with the catenary angle in the XGYG plane, given by the following expression:(5)ϕcat=arctan(x(N)-x(1)y(N)-y(1)),RZcat=[cos⁡(ϕcat)sin⁡(ϕcat)0-sin⁡(ϕcat)cos⁡(ϕcat)0001].The robot can move freely along the riser in theZl direction, but it was constrained in the other directions because it embraced the riser on the outside. If no water current was present, the nominal robot velocity propelled by a pair of thrusters should be ∥Vr∥ in the direction of Zl. However, due to the presence of the current, the absolute velocity of the robot was found by considering the sea current velocity component that is projected over the robot’s trajectory, which can change the robot’s progression velocity:(6)Vr_abs=-Vr+Vcur⋅Zl,
where Vr_abs is the absolute velocity and · the dot product. From a preliminary study, ∥Vr∥ was estimated as 1 m/s. Figure 6 illustrates robot Zl velocity in the particular sea current and thrusters conditions adopted in this paper.Figure 6
Velocity of the robot inZl direction considering a linear sea current profile.The effect of the longitudinal sheer between the robot and the riser due to the transversal current was not taken into account. This current component was expected essentially to increase the normal force that the robot applied to the outer surface of the riser. Because the riser was tightly fitted among the robot’s rigid free wells to avoid longitudinal and torsional slipping, the increase in the shear force of the wheels that could decelerate the robot was considered to be negligible.Since all the elements of the FEM mesh have approximately the same length, the time steps are no longer uniformly distributed with such variable velocity profile. The resulting variable time array was calculated by:(7)tvar(i)=‖D(i)‖‖Vr‖,i=1,…,N.This nonuniform time array was inconvenient for future calculations of velocity and acceleration profiles, and a new set of(x(i),y(i),z(i)) was found by spline interpolation using a uniform time array with the same limits of tvar. In the sequence, the local reference systems (Xl, Yl, Zl), the DCM matrix and the Euler angles were recalculated using this new set of coordinates, which were sampled by a uniform time array. All profiles were resampled at 5 Hz by spline interpolation, such that the high frequency riser motions could be followed (next subsection). The total mission time was Tm=1526s, and a total of M=Tm(s)×5(Hz)=7630 points were generated.


### 4.2. Riser Motion Effect

The dynamic FEM analysis allows finding, along the time, the altered geometry of the riser, with respect to the nominal static profile. The finite elements model had 170 beam elements with flexion, axial, and torsional deformations. Each element had a nominal length of 10 m. Dynamic FEM analysis provided theXG, YG, and ZG coordinates of only ten nodes with respect to time at fs=20Hz sampling rate. The analysis was run for a 100 second time span, but the first 50 seconds were disregarded to avoid the transient effect of the FEM solution. The second half of the time window was then replicated to reach the total mission time. Thus, the offset, phase, and amplitude parameters were preserved. This adjustment provided a matrix of ten lines (one for each node) and Tm×fs columns, which was resampled by successive spline interpolations. The first interpolation reduced the number of columns to M, the second expanded the lines up to the original FEM mesh (171 nodes) and the third resampled the lines again up to M. Therefore, three M×M perturbation matrices Perx, Pery, and Perz were obtained, one for each coordinate Xl, Yl, and Zl. In these matrices, each row was a particular riser deviation from the nominal profile and each column was a time step. The three components of a displacement vector d(i) were defined for each node i:(8)d(i)=[Per⁡x(i,i)Per⁡x(i,i)Per⁡x(i,i)],i=1,…,M.This vector was decomposed into its normaldn and tangential dt parts. Because the robot was free to move along the riser, only the normal component of d vector was effectively transmitted to the robot:(9)dt=(d⋅Du)Du,dn=d-dt.Therefore, the new coordinatesx, y, and z (in the local reference frame) of the riser path were given by the following sequence of Pi1¯new points, i=1,(10)Pi1¯new=Pi1¯nominal+dn.


### 4.3. Vortex-Induced Vibration (VIV) Effect

The sea current passing through a circular cylinder produces vortex-shedding in the wake, which causes the structure to vibrate. This complex fluid-structure interaction is called Vortex Induced Vibration (VIV), and it occurs predominantly on the cross-flow direction [16]. Simulating this effect is an arduous numerical problem. In this paper, we used experimental data obtained from a scale-model available at the Open-Source VIV Data Repository of the Center for Ocean Engineering at MIT [17]. Cross-flow displacement data that is available for tests performed in a bare cylinder, 20 mm diameter and 10 m length, which was donated by ExxonMobil, was used in our simulations. Several fluid velocities and both linearly sheared and linear flow conditions may be used. We chose a strong condition of regular flow, of approximately 1 m/s, which provided greater displacements compared to the sheared flow for the same nominal velocity.To adapt the experimental data to the riser that was being analyzed, the displacement was scaled by the test riser diameter and multiplied by the actual riser diameter. The frequency of shedding (fst) in cylinders with cross-flow is given by the following:(11)fst=StUD,
where St is the Strouhal number, U is the fluid velocity, and D is the cylinder diameter. Keeping St fixed, the ratio of fst between the experiment and the riser was 14.78. This factor was used to scale the time vector, such that the frequency vector of the data spectrum, which was used to simulate riser VIV, was divided by this quantity. The set of ten points along the cylinder where displacements were measured was associated with the closest nodes of the FEM model. A window of data, without transient effects, was replicated ten times until the total mission time was achieved, similar to the previous section, and the same interpolation procedure was applied. Finally, the displacement caused by VIV was rotated to become perpendicular to the catenary plane and added to (10).Figure7 shows the RMS profiles of the riser displacements in the XG, YG, and ZG directions as a function of the normalized riser length. The figure also shows the RMS of VIV perpendicular to the catenary plane. The number 0 is for the wheel head, and 1 corresponds to the turret.RMS of riser motion inXl, Yl, and Zl directions and VIV perpendicular to the catenary plane.
(a)(b)(c)(d)


### 4.4. Sensor Velocity and Acceleration Profiles

To simulate IMU output, the expected physical acceleration at the sensor installation point must be found. This acceleration was used to feed the sensor model to find realistic sensor signals. The acceleration at pointPS expressed in the global reference frame was given by the following well-known kinematic equation:(12)aPSG=[ẍÿz̈]+α×PS+ω×ω×PS.The local accelerationaPSl and Coriolis terms are zero because the sensor is fixed in the robot body. α is the angular acceleration and ω is the angular velocity. Angular velocity expressed in the global reference frame was calculated by the methods used in [14, 18]:(13)ω=RzRy[ψ̇00]+Rz[0θ̇0]+[00ϕ̇].To solve (12) and (13), the time derivatives of displacement and rotation ẋ, ẍ, ẏ, ψ̇, and so forth were found by finite differences and considering all variables in the global reference system.


## 4.1. Sea Current Effect

A linear current profileVcur from the maximum sea current velocity (Vcurmax) at sea level to 0 at seabed [15] was considered as being aligned with the riser catenary plane. Vcurmax was assumed as 1.68 m/s, from Campos Basin data:(4)Vcur(i)=RZcat[0Vcurmax(z(N-1)-z(i))z(N-1)-z(1)0]i=1,…,N,
where RZcat is the rotation matrix associated with the catenary angle in the XGYG plane, given by the following expression:(5)ϕcat=arctan(x(N)-x(1)y(N)-y(1)),RZcat=[cos⁡(ϕcat)sin⁡(ϕcat)0-sin⁡(ϕcat)cos⁡(ϕcat)0001].The robot can move freely along the riser in theZl direction, but it was constrained in the other directions because it embraced the riser on the outside. If no water current was present, the nominal robot velocity propelled by a pair of thrusters should be ∥Vr∥ in the direction of Zl. However, due to the presence of the current, the absolute velocity of the robot was found by considering the sea current velocity component that is projected over the robot’s trajectory, which can change the robot’s progression velocity:(6)Vr_abs=-Vr+Vcur⋅Zl,
where Vr_abs is the absolute velocity and · the dot product. From a preliminary study, ∥Vr∥ was estimated as 1 m/s. Figure 6 illustrates robot Zl velocity in the particular sea current and thrusters conditions adopted in this paper.Figure 6
Velocity of the robot inZl direction considering a linear sea current profile.The effect of the longitudinal sheer between the robot and the riser due to the transversal current was not taken into account. This current component was expected essentially to increase the normal force that the robot applied to the outer surface of the riser. Because the riser was tightly fitted among the robot’s rigid free wells to avoid longitudinal and torsional slipping, the increase in the shear force of the wheels that could decelerate the robot was considered to be negligible.Since all the elements of the FEM mesh have approximately the same length, the time steps are no longer uniformly distributed with such variable velocity profile. The resulting variable time array was calculated by:(7)tvar(i)=‖D(i)‖‖Vr‖,i=1,…,N.This nonuniform time array was inconvenient for future calculations of velocity and acceleration profiles, and a new set of(x(i),y(i),z(i)) was found by spline interpolation using a uniform time array with the same limits of tvar. In the sequence, the local reference systems (Xl, Yl, Zl), the DCM matrix and the Euler angles were recalculated using this new set of coordinates, which were sampled by a uniform time array. All profiles were resampled at 5 Hz by spline interpolation, such that the high frequency riser motions could be followed (next subsection). The total mission time was Tm=1526s, and a total of M=Tm(s)×5(Hz)=7630 points were generated.


## 4.2. Riser Motion Effect

The dynamic FEM analysis allows finding, along the time, the altered geometry of the riser, with respect to the nominal static profile. The finite elements model had 170 beam elements with flexion, axial, and torsional deformations. Each element had a nominal length of 10 m. Dynamic FEM analysis provided theXG, YG, and ZG coordinates of only ten nodes with respect to time at fs=20Hz sampling rate. The analysis was run for a 100 second time span, but the first 50 seconds were disregarded to avoid the transient effect of the FEM solution. The second half of the time window was then replicated to reach the total mission time. Thus, the offset, phase, and amplitude parameters were preserved. This adjustment provided a matrix of ten lines (one for each node) and Tm×fs columns, which was resampled by successive spline interpolations. The first interpolation reduced the number of columns to M, the second expanded the lines up to the original FEM mesh (171 nodes) and the third resampled the lines again up to M. Therefore, three M×M perturbation matrices Perx, Pery, and Perz were obtained, one for each coordinate Xl, Yl, and Zl. In these matrices, each row was a particular riser deviation from the nominal profile and each column was a time step. The three components of a displacement vector d(i) were defined for each node i:(8)d(i)=[Per⁡x(i,i)Per⁡x(i,i)Per⁡x(i,i)],i=1,…,M.This vector was decomposed into its normaldn and tangential dt parts. Because the robot was free to move along the riser, only the normal component of d vector was effectively transmitted to the robot:(9)dt=(d⋅Du)Du,dn=d-dt.Therefore, the new coordinatesx, y, and z (in the local reference frame) of the riser path were given by the following sequence of Pi1¯new points, i=1,(10)Pi1¯new=Pi1¯nominal+dn.


## 4.3. Vortex-Induced Vibration (VIV) Effect

The sea current passing through a circular cylinder produces vortex-shedding in the wake, which causes the structure to vibrate. This complex fluid-structure interaction is called Vortex Induced Vibration (VIV), and it occurs predominantly on the cross-flow direction [16]. Simulating this effect is an arduous numerical problem. In this paper, we used experimental data obtained from a scale-model available at the Open-Source VIV Data Repository of the Center for Ocean Engineering at MIT [17]. Cross-flow displacement data that is available for tests performed in a bare cylinder, 20 mm diameter and 10 m length, which was donated by ExxonMobil, was used in our simulations. Several fluid velocities and both linearly sheared and linear flow conditions may be used. We chose a strong condition of regular flow, of approximately 1 m/s, which provided greater displacements compared to the sheared flow for the same nominal velocity.To adapt the experimental data to the riser that was being analyzed, the displacement was scaled by the test riser diameter and multiplied by the actual riser diameter. The frequency of shedding (fst) in cylinders with cross-flow is given by the following:(11)fst=StUD,
where St is the Strouhal number, U is the fluid velocity, and D is the cylinder diameter. Keeping St fixed, the ratio of fst between the experiment and the riser was 14.78. This factor was used to scale the time vector, such that the frequency vector of the data spectrum, which was used to simulate riser VIV, was divided by this quantity. The set of ten points along the cylinder where displacements were measured was associated with the closest nodes of the FEM model. A window of data, without transient effects, was replicated ten times until the total mission time was achieved, similar to the previous section, and the same interpolation procedure was applied. Finally, the displacement caused by VIV was rotated to become perpendicular to the catenary plane and added to (10).Figure7 shows the RMS profiles of the riser displacements in the XG, YG, and ZG directions as a function of the normalized riser length. The figure also shows the RMS of VIV perpendicular to the catenary plane. The number 0 is for the wheel head, and 1 corresponds to the turret.RMS of riser motion inXl, Yl, and Zl directions and VIV perpendicular to the catenary plane.
(a)(b)(c)(d)


## 4.4. Sensor Velocity and Acceleration Profiles

To simulate IMU output, the expected physical acceleration at the sensor installation point must be found. This acceleration was used to feed the sensor model to find realistic sensor signals. The acceleration at pointPS expressed in the global reference frame was given by the following well-known kinematic equation:(12)aPSG=[ẍÿz̈]+α×PS+ω×ω×PS.The local accelerationaPSl and Coriolis terms are zero because the sensor is fixed in the robot body. α is the angular acceleration and ω is the angular velocity. Angular velocity expressed in the global reference frame was calculated by the methods used in [14, 18]:(13)ω=RzRy[ψ̇00]+Rz[0θ̇0]+[00ϕ̇].To solve (12) and (13), the time derivatives of displacement and rotation ẋ, ẍ, ẏ, ψ̇, and so forth were found by finite differences and considering all variables in the global reference system.


## 5.  Instrumentation Model

A low-cost microstrain 3DMGX-1 IMU was the project choice. This is a compact and integrated device, suitable for a high-depth submarine application, where the electronic case must be as slender as possible, for mechanical structural reasons. It delivers 3D accelerations, angular velocities, and attitude/orientation matrix in a single-serial channel. The error characteristics of each output were modeled as a wide-band noise plus a first order moving bias Markov process [19]:(14)Y=(1+SFE)U+OFS+b+wbw,
whereY: simulated corrupted sensor output,U: simulated clean physical signal,SFE: Scale Factor Error,OFS: Offset Error,b: moving bias, 1st order Markov process,wbw: wide-band sensor noise.The wide-band sensor noise was defined aswbw=σbwν, ν=N[0,1] (white noise, zero mean). σbw2 is the wide-band noise variance. Scale factor error, 0.5% by sensor specifications, was considered as follows:(15)SFE=0.005×sgn⁡(U(-1,1)),
such that U(-1,1) is a uniform random variable in [-1,1] and sgn is the signum function. This equation introduces a random error limited to 0.5% in the sensor output signal amplitude, when substituted in (14).The moving bias was found by integrating the following finite difference equation with the Euler method:(16)bn+1=bn+Δt((-1/τr)bn+σsensor_biasν),σsensor_bias=2fsσbias2τr,
whereΔt: sampling period (s),fs: sampling frequency (Hz),σsensor_bias2: bias input noise variance,σbias2: random walk variance,τr: time constant (s).


### 5.1. Determination of Sensor Parameters

An experiment was performed to collect the error characteristics by keeping the sensor stationary in laboratory while recording signals. Total acquisition time was 4 h 21 min, with a 75 Hz sampling frequency, after the internal temperature was stabilized. The first determined parameter was the time constantτr. Two techniques were employed: the Allan variance plot [20] and the autocorrelation function [19]. Both techniques gave different results (see Table 2). According to the likeness between the real and the simulated signals, one of the two time-constant techniques was selected. For the attitude sensor, τr determined by Allan variance was chosen, whereas for acceleration and angular velocity, autocorrelation provided the best results.Except for angular velocity (AngRate), the experiment time series was low-pass filtered before calculatingσbias, to better characterize the low-frequency component of the error. For σbw, the high-frequency part was used. The cut-off frequencies were chosen by trial and error, and the agreement between the simulated and original error signals was observed for each case. Time series, Allan variance plot [21] and fast Fourier transform (FFT) were used to test the agreement between original and simulated error signals. In the case of attitude simulation, the Euler angles were discretized in the interval comprising -180.0328° to 180°, with steps of 0.1°, corresponding to sensor resolution. For acceleration, the resolution of the sensor was 0.2 mg, and the angular velocity was 0.01°/s. However, amplitude discretization was not implemented. Figure 8 shows the results for Pitch angle and Figure 9 shows the results for acceleration in the Xl direction. For other angles, directions and angular velocities, the results were similar.Comparative results for real and simulated IMU sensor error for Pitch angle. (a) Window of the time series (degrees). (b) Allan variance plot (degrees). (c) FFT of time series (degrees2).
(a)(b)(c)Comparative results for real and simulated IMU sensor error forXl acceleration, in g (Earth acceleration of gravity at sea level). (a) Window of the time series (g). (b) Allan variance plot (g). (c) FFT of time series (g2).
(a)(b)(c)The water pressure-based depth sensor specified for this robot was the Digiquartz 8CB4000-I Depth Sensor (Paroscientific, Inc., Redmond, WA), provided a 0.01% accuracy and a resolution of10-8 m. In the simulations, we considered that the depth measures were corrupted by the addition of a Gaussian noise with σ = 0.2 meters, according to the estimates made by Jalving [22].


## 5.1. Determination of Sensor Parameters

An experiment was performed to collect the error characteristics by keeping the sensor stationary in laboratory while recording signals. Total acquisition time was 4 h 21 min, with a 75 Hz sampling frequency, after the internal temperature was stabilized. The first determined parameter was the time constantτr. Two techniques were employed: the Allan variance plot [20] and the autocorrelation function [19]. Both techniques gave different results (see Table 2). According to the likeness between the real and the simulated signals, one of the two time-constant techniques was selected. For the attitude sensor, τr determined by Allan variance was chosen, whereas for acceleration and angular velocity, autocorrelation provided the best results.Except for angular velocity (AngRate), the experiment time series was low-pass filtered before calculatingσbias, to better characterize the low-frequency component of the error. For σbw, the high-frequency part was used. The cut-off frequencies were chosen by trial and error, and the agreement between the simulated and original error signals was observed for each case. Time series, Allan variance plot [21] and fast Fourier transform (FFT) were used to test the agreement between original and simulated error signals. In the case of attitude simulation, the Euler angles were discretized in the interval comprising -180.0328° to 180°, with steps of 0.1°, corresponding to sensor resolution. For acceleration, the resolution of the sensor was 0.2 mg, and the angular velocity was 0.01°/s. However, amplitude discretization was not implemented. Figure 8 shows the results for Pitch angle and Figure 9 shows the results for acceleration in the Xl direction. For other angles, directions and angular velocities, the results were similar.Comparative results for real and simulated IMU sensor error for Pitch angle. (a) Window of the time series (degrees). (b) Allan variance plot (degrees). (c) FFT of time series (degrees2).
(a)(b)(c)Comparative results for real and simulated IMU sensor error forXl acceleration, in g (Earth acceleration of gravity at sea level). (a) Window of the time series (g). (b) Allan variance plot (g). (c) FFT of time series (g2).
(a)(b)(c)The water pressure-based depth sensor specified for this robot was the Digiquartz 8CB4000-I Depth Sensor (Paroscientific, Inc., Redmond, WA), provided a 0.01% accuracy and a resolution of10-8 m. In the simulations, we considered that the depth measures were corrupted by the addition of a Gaussian noise with σ = 0.2 meters, according to the estimates made by Jalving [22].


## 6.  Localization Algorithm

The localization problem consisted of estimating a set of states, which included robot position, by reading accelerations as inputs and considering the depth sensor signals as the external measurements. Since a strapdown navigation scheme is used, the localization algorithm assumes that the accelerations measured in the local frame are transformed into the global frame, by using the rotation matrices calculated from Euler angles measured by the IMU. It is also assumed that if a DC-acceleration sensible accelerometer is being used, the vertical acceleration was compensated for gravity. The problem state equations are formulated as [23, 24]:(17)Xk+1=AXk+BUk+Cwk,Zk=HXk+vk.The state vectors and matrices are given by the following:(18)Xk=[PxPyPzVxVyVzδaxδayδaz]kT,Uk=[axayaz]kT,wk=[w1w2w3w4w5w6w7w8w9]kT,vk=[ν1ν2ν3]kT.A[9×9] is an identity matrix, such that the terms (1,4), (2,5), and (3,6) are equal to T; (4,7), (5,8), and (6,9) assumes the value -T; (1,7), (2,8), and (3,9) elements are -1/2T2. B[9×3] is a matrix of zeros except for the terms (1,1), (2,2), and (3,3) that are equal to 1/2T2; (4,1), (5,2), and (6,3) elements are T.H[3×9] is a matrix of zeros, with the exception of the (i,i) terms equal to 1. In such expressions, k is the sample number, X is the state vector, U the input vector, w is the process noise, P is the position, V is the velocity, a is acceleration, δa is the acceleration error, T is sampling period, Z is the external measurements vector, v is the sensor noise, and σ2 is the variance of each associated variable. The observability matrix of this system had a full-rank, that is, 9.An estimateX̂k of the state vector is found by a Kalman filter.Prediction:(19)X̂k-=AX̂k-1-+BUk,Pk-=APk-1-AT+Q.Update:(20)Kk=Pk-HT(HPk-HT+R)-1,X̂k=X̂k-+Kk(Zk-HX̂k-),Pk=(I-KkH)Pk-,
where P is the state estimate error covariance matrix, R is the sensor noise covariance matrix, Q is the process noise matrix, and K is the Kalman gain. This is a standard, simple and generic implementation of a Kalman filter estimator. More sophisticated approaches could be tested as well, although the results obtained with the standard KF were satisfactory for the intended application, as will be addressed in Section 7. Unlike the kinematic models that are usually used for land robots [25] or those built to operate over flat surfaces [26], the state-space model used here to formulate the KF does not use any a priori information of the vehicle characteristics. Nor the dynamical characteristics of the sensors are included in this KF formulation. Thus, in principle, the localization algorithm proposed here could be applied for similar applications with different instrumentation, robot, and environmental characteristics. However, some tuning work on the covariance matrices should likely to be necessary.


### 6.1. Processing of the External Measurements for the KF

Here, we focus the analysis on the key aspect of this particular localization problem, which is determining the external measurements vectorZk=[OkxGOkyGOkzG]T. OkxG, and so forth, are the coordinates of the robot position in the global reference frame as measured by an independent absolute sensor. The sensor, that is, in practice, available to perform such measurements is the depth sensor, which can only measure the ZG coordinate. The attempts to use this measurement exclusively in the Update phase of the KF gave nondrifting estimates in the vertical direction, as expectable, but not in the horizontal one.We proposed the following method to estimate the completeZk vector using data from Zk-1, OkzG (depth sensor signal), and the Euler angles measured by the IMU. Considering Figure 10, the vector joining Zk to Zk-1, in the local reference frame centered in Zk, is given by [00L]T. This vector lies in Zl direction of the local reference frame 1k, with modulus L, and coinciding with the instantaneous riser path. In the global reference frame, Zk-1 can be found as:(21)R1kG=RzRyRx(ψ,θ,ϕ)=[R11R12R13R21R22R23R31R32R33],(22)Zk-1=Zk-R1kGZk-11k=Zk-[R13LR23LR33L]T.Figure 10
Position external measurementsZk-1 and Zk and the corresponding local reference frames (Xl, Yl, Zl). L is the scalar displacement between the two positions.Thus,(23)(Zk-Zk-1)=[R13LR23LR33L]T.The same difference calculated in (23) can be expressed by a vector of global coordinate differences:(24)(Zk-Zk-1)=[ΔXΔYΔZ]T.SinceR33 and ΔZ are known, by comparing (23) and (24), it is possible to find L as:(25)L=ΔZR33.SubstitutingL in (23), the KF external measurements vector Zk is found. Alternatively, the displacement L could be measured by an odometer. However, in this case, special constructive care should be taken into account to prevent slippage. In any case, this approach assumes a previous observation in the Correction phase of the KF. However, KF theory assumes conditional independence among the observations. Therefore, only suboptimal estimation performance is likely to be expected.


## 6.1. Processing of the External Measurements for the KF

Here, we focus the analysis on the key aspect of this particular localization problem, which is determining the external measurements vectorZk=[OkxGOkyGOkzG]T. OkxG, and so forth, are the coordinates of the robot position in the global reference frame as measured by an independent absolute sensor. The sensor, that is, in practice, available to perform such measurements is the depth sensor, which can only measure the ZG coordinate. The attempts to use this measurement exclusively in the Update phase of the KF gave nondrifting estimates in the vertical direction, as expectable, but not in the horizontal one.We proposed the following method to estimate the completeZk vector using data from Zk-1, OkzG (depth sensor signal), and the Euler angles measured by the IMU. Considering Figure 10, the vector joining Zk to Zk-1, in the local reference frame centered in Zk, is given by [00L]T. This vector lies in Zl direction of the local reference frame 1k, with modulus L, and coinciding with the instantaneous riser path. In the global reference frame, Zk-1 can be found as:(21)R1kG=RzRyRx(ψ,θ,ϕ)=[R11R12R13R21R22R23R31R32R33],(22)Zk-1=Zk-R1kGZk-11k=Zk-[R13LR23LR33L]T.Figure 10
Position external measurementsZk-1 and Zk and the corresponding local reference frames (Xl, Yl, Zl). L is the scalar displacement between the two positions.Thus,(23)(Zk-Zk-1)=[R13LR23LR33L]T.The same difference calculated in (23) can be expressed by a vector of global coordinate differences:(24)(Zk-Zk-1)=[ΔXΔYΔZ]T.SinceR33 and ΔZ are known, by comparing (23) and (24), it is possible to find L as:(25)L=ΔZR33.SubstitutingL in (23), the KF external measurements vector Zk is found. Alternatively, the displacement L could be measured by an odometer. However, in this case, special constructive care should be taken into account to prevent slippage. In any case, this approach assumes a previous observation in the Correction phase of the KF. However, KF theory assumes conditional independence among the observations. Therefore, only suboptimal estimation performance is likely to be expected.


## 7. Simulation Results

### 7.1. Static Results

In the static case no riser motion is present. Thestandard position was used as the basis to generate the dynamic analysis and is an intermediary neutral configuration between the near (ship displaces towards wheel head) and the far (the opposite direction) positions. Figure 11 shows the robot trajectory and attitude, represented by the successive local reference frames, for the neutral position. In Figure 12, the associated localization errors are shown, as a function of distance from the seabed (using a grid of 1.680 points, or 1.44 Hz). The absolute errors associated with the neutral, far, and near configurations are shown in Figure 13, in a semilog scale. To study the impact in the localization error of the grid refinement grade, an additional error curve was generated, for the neutral configuration, with a 7,630 points (5 Hz) grid, the same used for the dynamic analysis.Figure 11
Problem neutral configuration, showing the FPSO, the riser profile and the robot attitude. The orthogonal arrows express the successive local reference frames position and orientation along the riser path. Only 1/8 of the 171 FEM mesh nodes used in the static analysis is shown to facilitate visualization.Localization error (in global coordinates) for the static neutral case as a function of the distance from the seabed.
(a)(b)Total localization error (in global coordinates) as a function of distance from the seabed, for neutral, near and far configurations with the mesh of 1.680 points (1.44 Hz). The neutral configuration error for the mesh with 7.630 points (5 Hz) is also shown.
(a)(b)Figure12 shows that the greatest source of error occurs in the XG direction, while in the deepness coordinate ZG, the error is low. In addition, the total error increased as the robot tilted horizontally. This trend was expected due to the relative loss of accuracy of the depth sensor to estimate the L parameter when using (21)–(25). When the depth measure difference ΔZ between the two samples is small, a decrease in signal to Gaussian noise ratio is observed. From an application point of view, the robot should not work too close to the point where the riser periodically touches the seabed, the touch down point (TDP). This safety margin is established as being of 50 m. At this point, the expected positioning error for the static case was approximately 10 m. The error was slightly higher in the far case (Figure 13), where the riser was more horizontal than in the near condition. By increasing the sampling frequency (Figure 12), the error was likely to decrease in most of the path, but when the robot approximates to the TDP, the localization system fails.


### 7.2. Dynamic Results

By including the wave and VIV effects in the riser nodal displacement profile, the robot trajectory becomes more complex (Figure14), as expectable. The associated localization errors are shown in Figure 15, using the 5 Hz mesh. Figure 16 shows a zoom in the XGYG plane, close to TDP, of both the robot trajectory and the external measurements, calculated with (21)–(25), as well as the Kalman filter position states estimations. The KF estimated the trajectory accurately, up to the point where the processed external measurements lost accuracy, below 15 m distance to seabed. The KF estimates essentially followed the external measurements, in part because the KF Update frequency was the same as the Prediction. By comparing the errors in both static and dynamic cases (Figure 17), analogous profiles were found. It means that the Kalman filter was able to closely track the robot’s trajectory, even in the dynamic case, and the loss of algorithm localization performance was mainly caused, after all, by depth sensor differential signal to noise ratio.Figure 14
Robot attitude profile in the dynamic condition, including wave and VIV effects. Only 1/150 of the 7.630 grid points are shown to keep the figure clearer.Localization error (in global coordinates) as a function of the distance from the seabed, for the dynamic condition.
(a)(b)Figure 16
XGZG plane close to TDP for the dynamic condition, showing simulated robot trajectory (dashed), simulated and processed external measurements from depth sensor and IMU models (Δ) and Kalman Filter estimates (clearer dots).Figure 17
Mean localization error (in global coordinates) as a function of distance from the seabed, for both static and dynamic conditions.


## 7.1. Static Results

In the static case no riser motion is present. Thestandard position was used as the basis to generate the dynamic analysis and is an intermediary neutral configuration between the near (ship displaces towards wheel head) and the far (the opposite direction) positions. Figure 11 shows the robot trajectory and attitude, represented by the successive local reference frames, for the neutral position. In Figure 12, the associated localization errors are shown, as a function of distance from the seabed (using a grid of 1.680 points, or 1.44 Hz). The absolute errors associated with the neutral, far, and near configurations are shown in Figure 13, in a semilog scale. To study the impact in the localization error of the grid refinement grade, an additional error curve was generated, for the neutral configuration, with a 7,630 points (5 Hz) grid, the same used for the dynamic analysis.Figure 11
Problem neutral configuration, showing the FPSO, the riser profile and the robot attitude. The orthogonal arrows express the successive local reference frames position and orientation along the riser path. Only 1/8 of the 171 FEM mesh nodes used in the static analysis is shown to facilitate visualization.Localization error (in global coordinates) for the static neutral case as a function of the distance from the seabed.
(a)(b)Total localization error (in global coordinates) as a function of distance from the seabed, for neutral, near and far configurations with the mesh of 1.680 points (1.44 Hz). The neutral configuration error for the mesh with 7.630 points (5 Hz) is also shown.
(a)(b)Figure12 shows that the greatest source of error occurs in the XG direction, while in the deepness coordinate ZG, the error is low. In addition, the total error increased as the robot tilted horizontally. This trend was expected due to the relative loss of accuracy of the depth sensor to estimate the L parameter when using (21)–(25). When the depth measure difference ΔZ between the two samples is small, a decrease in signal to Gaussian noise ratio is observed. From an application point of view, the robot should not work too close to the point where the riser periodically touches the seabed, the touch down point (TDP). This safety margin is established as being of 50 m. At this point, the expected positioning error for the static case was approximately 10 m. The error was slightly higher in the far case (Figure 13), where the riser was more horizontal than in the near condition. By increasing the sampling frequency (Figure 12), the error was likely to decrease in most of the path, but when the robot approximates to the TDP, the localization system fails.


## 7.2. Dynamic Results

By including the wave and VIV effects in the riser nodal displacement profile, the robot trajectory becomes more complex (Figure14), as expectable. The associated localization errors are shown in Figure 15, using the 5 Hz mesh. Figure 16 shows a zoom in the XGYG plane, close to TDP, of both the robot trajectory and the external measurements, calculated with (21)–(25), as well as the Kalman filter position states estimations. The KF estimated the trajectory accurately, up to the point where the processed external measurements lost accuracy, below 15 m distance to seabed. The KF estimates essentially followed the external measurements, in part because the KF Update frequency was the same as the Prediction. By comparing the errors in both static and dynamic cases (Figure 17), analogous profiles were found. It means that the Kalman filter was able to closely track the robot’s trajectory, even in the dynamic case, and the loss of algorithm localization performance was mainly caused, after all, by depth sensor differential signal to noise ratio.Figure 14
Robot attitude profile in the dynamic condition, including wave and VIV effects. Only 1/150 of the 7.630 grid points are shown to keep the figure clearer.Localization error (in global coordinates) as a function of the distance from the seabed, for the dynamic condition.
(a)(b)Figure 16
XGZG plane close to TDP for the dynamic condition, showing simulated robot trajectory (dashed), simulated and processed external measurements from depth sensor and IMU models (Δ) and Kalman Filter estimates (clearer dots).Figure 17
Mean localization error (in global coordinates) as a function of distance from the seabed, for both static and dynamic conditions.


## 8.  Discussion and Conclusions

The simulations show that the robot trajectory using the localization algorithm, in the dynamic case, presented a wavy pattern. The pattern corresponds to the path that is traveled by the robot, and not the nominal riser profile, as expected. If an estimate of the riser catenary shape is desired, the obtained path could be used to fit a smooth profile curve.In the more vertical part of the riser (above 15 m from seabed), the average estimated error (standard dev.) was 0.76 (0.47) m (Figure15). This error can be considered sufficiently small for localizing the riser in the oil field and to feed nodal displacement constraints in a structural Finite Elements tension analysis. This accuracy is also satisfactory for localizing flaws detected by the NDTs.The localization algorithm depends strongly on the external measurements, which are provided by the depth sensor and the IMU, using (21)–(25). However, close to seabed, the simulated robot trajectory tilts towards a more horizontal attitude, and following the riser section which is closer to the TDP. In this part of the trajectory, Gaussian noise from the depth sensor more significantly corrupts the relative depth measurements. This effect increases when a period between two successive measurements becomes shorter. The 5 Hz is the minimum estimated rate necessary to reproduce VIV and the wave motion data. On the other hand, the minimum safety distance between the robot and the seabed coincides with the condition which the localization algorithm loose accuracy. As a possible extreme limit, at 15 m from the seabed, the expected error should be also around 15 m (Figure 15).In our opinion, the simple and standard KF implementation that was presented in this study is sufficient to provide the localization accuracy required for the proposed application. In the case of rigid risers and other predominantly vertical offshore structures, the observed loss of accuracy is not expected to occur. The kinematic model used to formulate the sensor fusion problem was linear, and therefore the standard KF implementation is appropriate. Other models could be proposed, incorporating robot, riser or sea dynamic features. In such cases, using an extended Kalman filter or a particle filter should be necessary for handling the associated nonlinearities. In any case, the external measurements reliability seems to be key aspect for localization accuracy in this problem.Using a higher grade low drift IMU instead of a low-cost one could eliminate the need for external measurements in a critical region but at the price of increasing payload, cost and volume. Acoustic positioning systems could be used to localize the robot in the most horizontal parts of the trajectory. However, these systems work at a very low sampling frequency (0.1 to 1 Hz) [8] and, at best, they could only probably cast the external motion envelope of the robot. More practically, improved accuracy could be achieved in the section close to the TDP by using an odometer to estimate L in (25). As such, the mechanical setup of an odometer should be accomplished very carefully to prevent slippage. Preinstalled RFID (Radio Frequency Identification) [27] tags or visual marks along the riser, signaling the actual length, are possibilities that could be explored. Switching the depth sensor with other kinds of sensors for external measurements, close to TDP, could be performed by using a fuzzy expert system [28]. In future studies, experimental laboratory and field tests should be performed, using an acoustic localization system to assess the bounds of the navigation errors estimated in this work.



---

*Source: 638394-2012-05-07.xml*