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.

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.

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, is contaminated 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 a severe 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.

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 Figure 3 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).

3. Architecture of the Localization System

The proposed localization system, that will be simulated numerically here, is shown in Figure 4. 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.

4. Generation of Static and Dynamic Sensor Profiles

A simulated profile of inertial sensor physical excitation was obtained. First, the nominal or static riser geometry was obtained from the FEM (finite element model) analysis. The global reference system (𝐗𝐆, 𝐘𝐆, 𝐙𝐆) was positioned in the turret center, such that 𝑋 axis points to an arbitrary direction (e.g., north), 𝑍 points up, and 𝑌 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.

The local reference system (𝐗𝐥, 𝐘𝐥, 𝐙𝐥) 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 𝑃𝐺𝑖=(𝑥(𝑖),𝑦(𝑖),𝑧(𝑖)), 𝑖=1,,𝑁 (𝑁 is number of nodes in the mesh and the superscript 𝐺 express the reference system). Each 𝐙𝐥 vector is found from a unit vector 𝐃𝐮 that connects node 𝑖 to node 𝑖+1 of the FEM mesh. The negative sign for defining 𝐙𝐥 in (2) is due to the fact that 𝐃𝐮 points down:,𝐃𝐃=𝑥(𝑖+1)𝑥(𝑖)𝑦(𝑖+1)𝑦(𝑖)𝑧(𝑖+1)𝑧(𝑖)𝐮=𝐃,𝐙𝐃(1)𝐥=𝐃𝐮.(2)

In the sequence, a particular 𝐘𝐥 was chosen (along with its orthogonal counterpart 𝐗𝐥), 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, 𝐘𝐥 is simultaneously orthogonal to a point sufficiently far in the 𝑋 direction and to −𝐙𝐥. The 𝐗𝐥 vector was orthogonal to both 𝐙𝐥 and 𝐘𝐥. For each local reference system 𝑖, a directors cosine matrix (DCM) was defined as follows:𝑋DCM=1𝑥𝑌1𝑥𝑍1𝑥𝑋1𝑦𝑌1𝑦𝑍1𝑧𝑋1𝑧𝑌1𝑧𝑍1𝑧,(3) where 𝐗𝐥=(𝑋1𝑥,𝑋1𝑦,𝑋1𝑧)𝑇 is the unit vector that defines 𝐗𝐥, and similarly 𝐘𝐥 and 𝐙𝐥.

Euler angles 𝜓, 𝜃, and 𝜙 were defined between the 𝐗𝐆 and 𝐗𝐥, 𝐘𝐆 and 𝐘𝐥, 𝐙𝐆 and 𝐙𝐥 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 𝑅𝐺1=𝑅𝑥𝑅𝑦𝑅𝑧 was posteriorly verified to be the same as DCM, in order to determine if singularities were present. A fixed point 𝑃𝑆=(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 𝑌 distance from the riser section center, where the origin of the local reference frame is positioned. The coordinates of 𝑃𝑆, expressed in the global coordinate frame, were assumed as the sensor displacement profile.

4.1. Sea Current Effect

A linear current profile 𝐕𝐜𝐮𝐫 from the maximum sea current velocity (𝑉curmax) at sea level to 0 at seabed [15] was considered as being aligned with the riser catenary plane. 𝑉curmax was assumed as 1.68 m/s, from Campos Basin data:𝐕𝐜𝐮𝐫(𝐢)=𝑅𝑍cat0𝑉curmax(𝑧(𝑁1)𝑧(𝑖))0𝑧(𝑁1)𝑧(1)𝑖=1,,𝑁,(4) where 𝑅𝑍cat is the rotation matrix associated with the catenary angle in the 𝑋𝐺𝑌𝐺 plane, given by the following expression:𝜙cat=arctan𝑥(𝑁)𝑥(1),𝑦(𝑁)𝑦(1)𝑅𝑍cat=𝜙coscat𝜙sincat0𝜙sincat𝜙coscat0.001(5)

The robot can move freely along the riser in the 𝐙𝐥 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 𝐕𝐫 in the direction of 𝐙𝐥. 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) where 𝐕𝐫_𝐚𝐛𝐬 is the absolute velocity and the dot product. From a preliminary study, 𝐕𝐫 was estimated as 1 m/s. Figure 6 illustrates robot 𝐙𝐥 velocity in the particular sea current and thrusters conditions adopted in this paper.

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:𝑡var(𝑖)=𝐃(𝐢)𝐕𝐫,𝑖=1,,𝑁.(7)

This nonuniform time array was inconvenient for future calculations of velocity and acceleration profiles, and a new set of (𝑥(𝑖),𝑦(𝑖),𝑧(𝑖)) was found by spline interpolation using a uniform time array with the same limits of 𝑡var. In the sequence, the local reference systems (𝐗𝐥, 𝐘𝐥, 𝐙𝐥), 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 𝑇𝑚=1526s, and a total of 𝑀=𝑇𝑚(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 the 𝐗𝐆, 𝐘𝐆, and 𝐙𝐆 coordinates of only ten nodes with respect to time at 𝑓𝑠=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 𝑇𝑚×𝑓𝑠 columns, which was resampled by successive spline interpolations. The first interpolation reduced the number of columns to 𝑀, the second expanded the lines up to the original FEM mesh (171 nodes) and the third resampled the lines again up to 𝑀. Therefore, three 𝑀×𝑀 perturbation matrices Per𝑥, Per𝑦, and Per𝑧 were obtained, one for each coordinate 𝐗𝐥, 𝐘𝐥, and 𝐙𝐥. 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 𝐝(𝐢) were defined for each node 𝑖:𝐝(𝐢)=Per𝑥(𝑖,𝑖)Per𝑥(𝑖,𝑖)Per𝑥(𝑖,𝑖),𝑖=1,,𝑀.(8)

This vector was decomposed into its normal 𝐝𝐧 and tangential 𝐝𝐭 parts. Because the robot was free to move along the riser, only the normal component of 𝐝 vector was effectively transmitted to the robot:𝐝𝐭𝐝=(𝐝𝐃𝐮)𝐃𝐮,𝐧=𝐝𝐝𝐭.(9)

Therefore, the new coordinates 𝑥, 𝑦, and 𝑧 (in the local reference frame) of the riser path were given by the following sequence of 𝑃1𝑖new points, 𝑖=1,𝑃1𝑖new=𝑃1𝑖nominal+𝐝𝐧.(10)

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 (𝑓st) in cylinders with cross-flow is given by the following:𝑓st=St𝑈𝐷,(11) where St is the Strouhal number, 𝑈 is the fluid velocity, and 𝐷 is the cylinder diameter. Keeping St fixed, the ratio of 𝑓st 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).

Figure 7 shows the RMS profiles of the riser displacements in the 𝐗𝐆, 𝐘𝐆, and 𝐙𝐆 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.

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 point 𝑃𝑆 expressed in the global reference frame was given by the following well-known kinematic equation:𝐚𝐏𝐒𝐆=̈𝑥̈𝑦̈𝑧+𝜶×𝑃𝑆+𝝎×𝝎×𝑃𝑆.(12)

The local acceleration 𝐚𝐏𝐒𝐥 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]:𝝎=𝑅𝑧𝑅𝑦00̇𝜓+𝑅𝑧0̇𝜃0+00̇𝜙.(13)

To solve (12) and (13), the time derivatives of displacement and rotation ̇𝑥, ̈𝑥, ̇𝑦, ̇𝜓, 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]:𝑌=(1+SFE)𝑈+OFS+𝑏+𝑤bw,(14) where𝑌: simulated corrupted sensor output,𝑈: simulated clean physical signal,SFE: Scale Factor Error,OFS: Offset Error,𝑏: moving bias, 1st order Markov process,𝑤bw: wide-band sensor noise.

The wide-band sensor noise was defined as 𝑤bw=𝜎bw𝜈, 𝜈=𝑁[0,1] (white noise, zero mean). 𝜎2bw is the wide-band noise variance. Scale factor error, 0.5% by sensor specifications, was considered as follows:SFE=0.005×sgn(𝑈(1,1)),(15) such that 𝑈(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:𝑏𝑛+1=𝑏𝑛+Δ𝑡1/𝜏𝑟𝑏𝑛+𝜎sensor_bias𝜈,𝜎sensor_bias=2𝑓𝑠𝜎2bias𝜏𝑟,(16) whereΔ𝑡: sampling period (s),𝑓𝑠: sampling frequency (Hz),𝜎2sensor_bias: bias input noise variance,𝜎2bias: random walk variance,𝜏𝑟: 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 𝜏𝑟. 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, 𝜏𝑟 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 𝐗𝐥 direction. For other angles, directions and angular velocities, the results were similar.

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 of 108 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]:𝐗𝑘+1=𝐴𝐗𝑘+𝐵𝐔𝑘+𝐶𝐰𝑘,𝐙𝑘=𝐻𝐗𝑘+𝐯𝑘.(17)

The state vectors and matrices are given by the following:𝐗𝑘=𝑃𝑥𝑃𝑦𝑃𝑧𝑉𝑥𝑉𝑦𝑉𝑧𝛿𝑎𝑥𝛿𝑎𝑦𝛿𝑎𝑧𝑇𝑘,𝐔𝑘=𝑎𝑥𝑎𝑦𝑎𝑧𝑇𝑘,𝐰𝑘=𝑤1𝑤2𝑤3𝑤4𝑤5𝑤6𝑤7𝑤8𝑤9𝑇𝑘,𝐯𝑘=𝜈1𝜈2𝜈3𝑇𝑘.(18)

𝐴[9×9] is an identity matrix, such that the terms (1,4), (2,5), and (3,6) are equal to 𝑇; (4,7), (5,8), and (6,9) assumes the value 𝑇; (1,7), (2,8), and (3,9) elements are 1/2𝑇2. 𝐵[9×3] is a matrix of zeros except for the terms (1,1), (2,2), and (3,3) that are equal to 1/2𝑇2; (4,1), (5,2), and (6,3) elements are 𝑇.

𝐻[3×9] is a matrix of zeros, with the exception of the (i,i) terms equal to 1. In such expressions, 𝑘 is the sample number, 𝐗 is the state vector, 𝐔 the input vector, 𝐰 is the process noise, 𝑃 is the position, 𝑉 is the velocity, 𝑎 is acceleration, 𝛿𝑎 is the acceleration error, 𝑇 is sampling period, 𝐙 is the external measurements vector, 𝐯 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 estimate 𝐗𝑘 of the state vector is found by a Kalman filter.

Prediction:𝐗𝑘𝐗=𝐴𝑘1+𝐵𝐔𝑘,𝑃𝑘=𝐴𝑃𝑘1𝐴𝑇+𝑄.(19)

Update:𝐾𝑘=𝑃𝑘𝐻𝑇(𝐻𝑃𝑘𝐻𝑇+𝑅)1,𝐗𝑘=𝐗𝑘+𝐾𝑘𝐙𝑘𝐗𝐻𝑘,𝑃𝑘=𝐼𝐾𝑘𝐻𝑃𝑘,(20) where 𝑃 is the state estimate error covariance matrix, 𝑅 is the sensor noise covariance matrix, 𝑄 is the process noise matrix, and 𝐾 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 vector 𝐙𝑘=[𝑂𝐺𝑘𝑥𝑂𝐺𝑘𝑦𝑂𝐺𝑘𝑧]𝑇. 𝑂𝐺𝑘𝑥, 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 𝐙𝐆 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 complete 𝐙𝑘 vector using data from 𝐙𝑘1, 𝑂𝐺𝑘𝑧 (depth sensor signal), and the Euler angles measured by the IMU. Considering Figure 10, the vector joining 𝐙𝑘 to 𝐙𝑘1, in the local reference frame centered in 𝐙𝑘, is given by [00𝐿]𝑇. This vector lies in 𝐙𝐥 direction of the local reference frame 1𝑘, with modulus 𝐿, and coinciding with the instantaneous riser path. In the global reference frame, 𝐙𝑘1 can be found as:𝑅𝐺1𝑘=𝑅𝑧𝑅𝑦𝑅𝑥𝑅(𝜓,𝜃,𝜙)=11𝑅12𝑅13𝑅21𝑅22𝑅23𝑅31𝑅32𝑅33𝐙,(21)𝑘1=𝐙𝑘𝑅𝐺1𝑘𝐙1𝑘𝑘1=𝐙𝑘𝑅13𝐿𝑅23𝐿𝑅33𝐿𝑇.(22)

Thus,𝐙𝑘𝐙𝑘1=𝑅13𝐿𝑅23𝐿𝑅33𝐿𝑇.(23)

The same difference calculated in (23) can be expressed by a vector of global coordinate differences:𝐙𝑘𝐙𝑘1=Δ𝑋Δ𝑌Δ𝑍𝑇.(24)

Since 𝑅33 and Δ𝑍 are known, by comparing (23) and (24), it is possible to find 𝐿 as:𝐿=Δ𝑍𝑅33.(25)

Substituting 𝐿 in (23), the KF external measurements vector 𝐙𝑘 is found. Alternatively, the displacement 𝐿 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. The standard 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 12 shows that the greatest source of error occurs in the 𝐗𝐆 direction, while in the deepness coordinate 𝐙𝐆, 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 𝐿 parameter when using (21)–(25). When the depth measure difference Δ𝑍 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 (Figure 14), as expectable. The associated localization errors are shown in Figure 15, using the 5 Hz mesh. Figure 16 shows a zoom in the 𝑋𝐺𝑌𝐺 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.

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 (Figure 15). 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 𝐿 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.

Acknowledgments

The authors are gratefully thankful to CAPES (Coordenacao de Aperfeicoamento de Pessoal de Nivel Superior), FINEP (Financiadora de Estudos e Projetos), CNPq (Conselho Nacional de Desenvolvimento Cientifico e Tecnologico), and FAPERJ (Fundacao de Amparo a Pesquisa do Estado do Rio de Janeiro) for financial support.