# Integration of a strapdown gravimeter system in an autonomous underwater vehicle

Clément Roussel, Jérôme Verdun, José Cali, Marcia Maia, Jean-François F d'Eu

► To cite this version:

Clément Roussel, Jérôme Verdun, José Cali, Marcia Maia, Jean-François F d'Eu. Integration of a strapdown gravimeter system in an autonomous underwater vehicle. ISPRS International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2015, XL-5/W5, pp.199 - 206. 10.5194/isprsarchives-xl-5-w5-199-2015. hal-03559412

HAL Id: hal-03559412

https://hal.science/hal-03559412

Submitted on 6 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# INTEGRATION OF A STRAPDOWN GRAVIMETER SYSTEM IN AN AUTONOMOUS UNDERWATER VEHICLE

C. Roussel<sup>a</sup>, J. Verdun<sup>a</sup>, J. Cali<sup>a</sup>, M. Maia<sup>b</sup>, J.F. d'EU<sup>b</sup>

## Commission V

**KEY WORDS:** Autonomous Underwater Vehicle, moving vectorial gravimetry and gradiometry, electrostatic accelerometer, Unscented Kalman filtering

## ABSTRACT:

We present a new mobile instrument for measuring dynamically the gravity vector and its gradients in underwater environment, named GRAVIMOB. Our instrument is a strapdown sensor, consisted of electrostatic accelerometers installed in a waterproof sphere. It is designed to be embedded in an Autonomous Underwater Vehicle (AUV). Since the positioning of an AUV is approximate in underwater environment, the key issue raised here is to estimate the uncertainty in the gravity field resulting from the use of such position data. This paper focuses on the assessment of the system performances. The latter have been studied by simulation with reference data calculated from actual submarine geological structures, on which different noise models have been added. Results show that spatial evolutions of the gravity field and statistical properties of stochastic processes affecting the measurements have to be considered carefully in order to minimize the error. The Unscented Kalman Filter (UKF) has been favored to the Extended Kalman Filter (EKF) by its ease of implementation and its better robustness to non-linearities.

## 1. INTRODUCTION

One of the major challenges of modern gravimetry consists in determining Earth's gravity field models covering the whole of wavelength range with the same reliability at local, regional and global spatial scales. Since the Earth's surface is covered by more than 70% by the oceans, gravity field measurement techniques must be adapted to the marine and submarine areas. Among these techniques, satellite altimetry enables to reach long wavelengths (≥ 100 km) and covers almost the entire surface of the oceans, but becomes inoperative at the land-sea interface where local effects such as ocean tides have a high amplitude. Airborne gravimetry, using airplanes or helicopters, or surface gravimetry using ships, permit the recovery of intermediate wavelengths (≥ 10 km) but requires the use of stabilized platforms because of the large velocity and attitude variations encountered with such carriers. Static measurements on the sea or ocean floor can be carried out very close to the geological structures, thus allowing short wavelengths (≤ 1 km) to be covered. However, they require the use of submersibles which must land on the bottom and remain stationary during the measuring time. Moreover, they cannot access to hilly zones and are very time-consuming.

A more recent technique able to measure the undersea gravity field consists in the integration of a measuring instrument in a Autonomous Underwater Vehicle (AUV). The main advantage of a mobile gravimetric measurement system operating in underwater environment is that the measure of the gravity field can be conducted very close to the geological sources along profiles of different depths where permanent instrumentation is not feasible due to the high seismic activity and volcanism such as encountered in subduction zones and oceanic ridges. Recovering the gravity field at different depths is very useful for interpreting gravity variations to infer the location and the shape of gravity sources. Consequently, many applications are possible in geodesy as well as in geophysics. Indeed, some marine geophysicists argue that the hydrothermal sources located in oceanic ridges should hold considerable mineral resources. Several studies have already been

carried out on the subject as in (Araya et al., 2011; Yan et al., 2012; Liu et al., 2010; Yao and Xiaorong, 2012). These studies focus on the vertical component of the gravity field and therefore require the use of stabilized platforms. But at a time when reduced size and limited available power of the submersibles become weighty arguments, this technique need to be rethought to make it compatible with current constraints.

We present here a new mobile instrument for measuring dynamically the gravity vector and its gradients in underwater environment, named GRAVIMOB. Our research relies on the experience gained in the field of land mobile gravimetry with the mobile gravimetry system LiMo-g (Light Moving gravimetry system) initiated by the Geodesy and Geomatics Laboratory (L2G) of the National Conservatory of Arts and Crafts (CNAM) and the Geodesy Laboratory (LAREG) of the French Geographic and Forest Institute (IGN) and studied during the doctoral thesis of Bertrand de Saint-Jean (2008). Our instrument is a strapdown sensor, composed of six electrostatic accelerometers divided into two triads, both installed in a waterproof sphere. It allows the three components of the gravity vector to be measured and does not require the use of a stabilized platform.

Obviously, the measurements performed with such an instrument only make sense if the movement of the AUV, its position and attitude, is known at all times. One of the main difficulties encountered in underwater environment is that the positioning of the carrier is approximate. The position, velocity, acceleration and attitude, essential parameters for the restitution of the submarine gravity field from accelerometer measurements, are obtained by combining data from different sensors – inertial navigation system and Doppler log – with a target position. This position is deduced from that of the accompanying ship, obtained by GNSS and transmitted to the submarine via acoustic signal. The key issue raised here is to estimate the uncertainty in the gravity field resulting from the use of such position data.

This paper focuses on the assessment of the system performances.

<footnote><sup>a</sup> Geodesy and Geomatics Laboratory (L2G)</footnote>

<footnote><sup>b</sup> Oceanic Domains Laboratory (LDO)</footnote>
Such performances have been evaluated thanks to numerical sim- ulations with synthetic data calculated from actual submarine ge- ological structures and actual AUV trajectories, obtained by ad- dition of different noises to reference data.

The paper is organized as follows: part two describes the instru- ment and the carrier in which it will be installed. Part three in- troduces the mathematical equations related to our instrument. Based on these equations, the four and the five parts expose the principle of the numerical simulation used for assessing the per- formances of our system. Reference data and noise models are also detailed. The results of various simulations will be given in part five. Part six introduces and justifies the choice of the Unscented Kalman Filter (UKF) as the processing method. Con- clusions and perspectives will be exhibited in part seven.

## 2. DESIGN

### 2.1 Sensors

Our instrumentation consists of six electrostatic accelerometers Q-FLEX QA-3000-020 developed by the Honeywell - AlliedSig- nal company (Figure 1). These six accelerometers are divided into two triads, named α and β, both installed in a waterproof glass sphere. Each triad allows the measurement of the three components of the specific force related to the points Mα and Mβ, both geometrically defined as being the respective intersec- tion point of the three sensitive axes of each triad. Note here that the frames defined by each triad are supposed orthogonal and co- incide with the frame of the carrier, called b-frame (Figure 2).

<figure><img src="image_2.png" /><figcaption>Figure 1: View of the Q-FLEX QA-3000-020 electrostatic ac- celerometer (source: asc-sensors.de)</figcaption></figure>

Technically, the operation of these sensors consists on keeping the equilibrium of a test mass mounted on a flexible blade made of quartz and maintained horizontal through the restoring force induced by the magnetic field generated by an electromagnet. The balance is controlled by a capacitive position sensor detect- ing any change of blade position. Under the effect of an acceler- ation, the latter is disturbed as the proof mass moves with a delay with respect to the housing of the sensor, due to its inertia and the bending of the blade. The capacitive position sensor detects the deviation between the proof mass and the frame and provides the electromagnet with the current necessary to bring back the flexi- ble blade to its equilibrium position. The force necessary for this compensation corresponds to the acceleration component under- gone in the sensitive direction of the accelerometer.

### 2.2 Carrier

GRAVIMOB is intended to be embedded in an Autonomous Un- derwater Vehicle (AUV) of the French Research Institute for Ex- ploitation of the Sea (IFREMER), named AsterX (Figure 3). This AUV is able to dive down to 3,000 meter depth and travel up to 100 km. Its total mass is 800 kg and the maximum weight of the scientific payload is 200 kg.

To ensure an autonomous navigation, this AUV is equipped with an Inertial Navigation System (INS) coupled to a Doppler Veloc-

<figure><img src="image_3.png" /><figcaption>Figure 2: Position of the two triads (α and β). The point P on this figure belongs to (Mα Mβ) and corresponds to the vehicle reference point. The sphere diameter is about 40 cm.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: The Autonomous Underwater Vehicle AsterX of the French Research Institute for Exploitation of the Sea (IFRE- MER). (source: http://flotte.ifremer.fr)</figcaption></figure>

ity Log (DVL). A surface ship, geolocated using Global Naviga- tion Satellite System (GNSS), monitors the submersible thanks to an Ultra Short Base Line (USBL) acoustic positioning system.

## 3. EQUATION OF MOVING-BASE GRAVIMETRY

### 3.1 Application of Newton’s Second Law

Let us consider one 3D accelerometer labelled by α. Let X̂i_α be the position vector of the proof mass inside the accelerome- ter triad expressed in the inertial frame abbreviated i-frame. The proof mass is located at the point Mα where acceleration mea- surement is performed. According to Newton’s Second Law ap- plied to the proof mass, the second-order time derivative X̂i_α of X̂i_α may be expressed as

$$ \ddot{X}_{\alpha}^{i} = g_{\alpha}^{i} + a_{\alpha}^{i} \qquad (1) $$

where ĝi_α is the gravitational acceleration and âi_α is the restoring force per unit of mass, that is the specific force exerted on the proof mass inside the sensor, both projected onto i-frame axes. The position of the point Mα is not directly accessible. Only the point P, previously defined (Figure 2), is assumed to be known. The position vector X̂i_α can therefore be written as:

$$ X_{\alpha}^{i} = C_{e}^{i} X_{P}^{e} + C_{b}^{i} L_{\alpha}^{b} \qquad (2) $$

where L̂b_α = PMb_α is the level arm between the point P located at the vehicle reference point and the measuring point Mα of the 3D accelerometer α, Ĉe_α is the rotation matrix between the earth frame, abbreviated e-frame and the i-frame and Ĉi_α is the rotation matrix between the b-frame and the i-frame. It must be empha- sized that the restoring force maintains the proof mass fixed with
respect to the case of the sensor which is itself fixed with respect to the b-frame. As a consequence, the components of the one column matrix $L_\alpha^b$ are constant independent of time. $C_e^i$ depends only on the speed of rotation of the Earth, $\omega_e$. $C_b^i$ may be decomposed as the product of the matrix $C_b^e$, the matrix $C_n^e$ and the matrix $C_e^i$:

$$C_b^i = C_e^i C_n^e C_b^n \quad (3)$$

where $C_n^e$ is the rotation matrix between the navigation frame, abbreviated n-frame and the e-frame depending only on the geographical coordinates of the point P (longitude $\lambda_P$, and latitude $\varphi_P$). $C_b^i$ is the rotation matrix between the b-frame and the n-frame depending only on the attitude of the carrier (heading, $\delta$, pitch, $\chi$ and roll, $\eta$).

The second-order time derivative of equation (2), combined with equation (1) leads to the basic equation of moving-base gravimetry:

$$g_{\alpha}^{n} = C_{e}^{n} \left[ \ddot{X}_{P}^{e} + 2 \Omega_{ie}^{e} \dot{X}_{P}^{e} + (\Omega_{ie}^{e} \Omega_{ie}^{e} + \dot{\Omega}_{ie}^{e}) X_{P}^{e} \right] + C_{b}^{n} (\Omega_{ib}^{b} \Omega_{ib}^{b} + \dot{\Omega}_{ib}^{b}) L_{\alpha}^{b} - C_{b}^{n} a_{\alpha}^{b} \quad (4)$$

where $\Omega_{ib}^b$ (resp. $\Omega_{ie}^e$) is the skew symmetric matrix associated with the rotation of the b-frame (resp. e-frame) with respect to the i-frame and $\dot{\Omega}_{ib}^b$ (resp. $\dot{\Omega}_{ie}^e$) its first-order time derivative, both expressed in the b-frame (resp. e-frame).

## 3.2 Simplifying assumptions

In the context of this article, we will ignore the influence of the rotation of the Earth with respect to the inertial frame. In other words, we consider that the inertial and the earth frames are the same. Therefore, the index i is replaced by the index e and equation $\Omega_{ie}^e = \dot{\Omega}_{ie}^e = 0$ holds. Equation (4) then becomes:

$$g_{\alpha}^{n} = C_{e}^{n} \ddot{X}_{P}^{e} + C_{b}^{n} (\Omega_{eb}^{b} \Omega_{eb}^{b} + \dot{\Omega}_{eb}^{b}) L_{\alpha}^{b} - C_{b}^{n} a_{\alpha}^{b} \quad (5)$$

Equation (5) can also be written at the point $M_\beta$ of the 3D acceleometer $\beta$ where acceleration measurement is also performed, thus giving:

$$g_{\beta}^{n} = C_{e}^{n} \ddot{X}_{P}^{e} + C_{b}^{n} (\Omega_{eb}^{b} \Omega_{eb}^{b} + \dot{\Omega}_{eb}^{b}) L_{\beta}^{b} - C_{b}^{n} a_{\beta}^{b} \quad (6)$$

Equations (5) and (6) show similarities. The acceleration term relating to the point P, and the factor multiplying the lever arms $L_\alpha^b$ and $L_\beta^b$ are the same in both cases. By linear combination of these two equations and under certain assumptions, it is possible to remove them.

## 3.3 Linear combinations

Gravity at $M_\alpha$ and $M_\beta$, can be related to the gravity at point P. Indeed, a first order approximation gives:

$$\left\{ \begin{array}{l} g_{\alpha}^{n} = g_{P}^{n} + \nabla g_{P}^{n} L_{\alpha}^{n} \\ g_{\beta}^{n} = g_{P}^{n} + \nabla g_{P}^{n} L_{\beta}^{n} \end{array} \right. , \quad (7)$$

where $\nabla g_P^n$ is the gravity gradient tensor at point P expressed in the n-frame and corresponding there to the 3 × 3 matrix:

$$\Delta g_{P}^{n} = \begin{bmatrix} \partial_{E} g^{E} & \partial_{N} g^{E} & \partial_{U} g^{E} \\ \partial_{E} g^{N} & \partial_{N} g^{N} & \partial_{U} g^{N} \\ \partial_{E} g^{U} & \partial_{N} g^{U} & \partial_{U} g^{U} \end{bmatrix} \quad (8)$$

where $g^E$, $g^N$ and $g^U$ are the 3 components of the gravity vector $g_P$ when expressed in the n-frame and by denoting $x_E$, $x_N$, $x_U$

the 3 coordinates of the point P in the n-frame, we have:

$$\partial_i g^j = \frac{\partial g^j}{\partial x_i} \quad (9)$$

for $i, j = E, N, U$.

Under the assumption that $L_\alpha = -L_\beta$ or in other words, that point P is the midpoint of the segment $[M_\alpha M_\beta]$, the sum of the previous two equations (7) leads to:

$$\frac{g_{\alpha}^{n} + g_{\beta}^{n}}{2} = g_{P}^{n} \quad (10)$$

By substituting $g_\alpha^n$ and $g_\beta^n$ using the expressions established in (5) and (6), we finally obtain:

$$g_{P}^{n} = C_{e}^{n} \ddot{X}_{P}^{e} - C_{b}^{n} \left( \frac{a_{\alpha}^{b} + a_{\beta}^{b}}{2} \right) \quad (11)$$

The same reasoning can be emitted by performing the difference of the equations listed in (7), we can then deduce:

$$\nabla g_{P}^{n} u_{\beta\alpha}^{n} = \frac{g_{\alpha}^{n} - g_{\beta}^{n}}{\| M_{\beta} M_{\alpha} \|} \quad (12)$$

and, according to (5) and (6):

$$\nabla g_{P}^{n} u_{\beta\alpha}^{n} = \frac{\left( C_{b}^{n} (\Omega_{eb}^{b} \Omega_{eb}^{b} + \dot{\Omega}_{eb}^{b}) M_{\beta} M_{\alpha}^{b} - C_{b}^{n} (a_{\alpha}^{b} - a_{\beta}^{b}) \right)}{\| M_{\beta} M_{\alpha} \|} \quad (13)$$

where $M_\beta M_\alpha$ is the vector connecting the points $M_\beta$ et $M_\alpha$, $\|M_\beta M_\alpha\|$ its norm and $u_{\beta\alpha}$ is unit vector defined as:

$$u_{\beta\alpha} = \frac{M_{\beta} M_{\alpha}}{\| M_{\beta} M_{\alpha} \|} \quad (14)$$

$\nabla g_P^n u_{\beta\alpha}^n$ can be viewed as the 3-component gravity gradient vector in the direction defined by the vector $u_{\beta\alpha}$.

# 4. NUMERICAL SIMULATION

## 4.1 Principle and purpose

In a simplified manner, equations (11) and (13) may be writing as follows:

$$f: \begin{pmatrix} \lambda_P, \varphi_P, h_P \\ \delta, \chi, \eta \\ a_\alpha, a_\beta \end{pmatrix} \rightarrow \left\{ \begin{array}{l} g_P \\ \nabla g_P u_{\beta\alpha} \end{array} \right. , \quad (15)$$

where $f$ is a multivariate function mapping $\mathbb{R}^{12}$ into $\mathbb{R}^6$, which relates the geodetic coordinates of the point P (longitude $\lambda_P$, latitude $\varphi_P$, ellipsoidal height $h_P$), the attitude angles (heading $\delta$, pitch $\chi$, roll $\eta$) and the 2 × 3 = 6 components of acceleration $a_\alpha$ and $a_\beta$ to the 3 components of the gravity vector $g_P$ and the 3 gravity gradients in the direction of vector $u_{\beta\alpha}$. Since the input data (position, attitude and restoring accelerations) are subject to uncertainties, bias, drift and other faults, we attempt here to determine the uncertainty affecting the output data, namely the gravity vector $g_P$ and the gravity gradient vector $\nabla g_P u_{\beta\alpha}$. As the function $f$ is highly non-linear, we opted for an assessment of the uncertainty of the output error by means of Monte Carlo statistical method. For this task, let us consider arbitrary noise models representing faults on the different inputs given by equation (15). For the i-th simulation, the new input variables can be
written as:

$$
\begin{align}
\tilde{\lambda}_{P,i} &= \lambda_P + \epsilon_{\lambda,i} \nonumber \\
\tilde{\varphi}_{P,i} &= \varphi_P + \epsilon_{\varphi,i} \nonumber \\
\tilde{h}_{P,i} &= h_P + \epsilon_{h,i} \nonumber \\
\tilde{\delta}_i &= \delta + \epsilon_{\delta,i} \tag{16} \\
\tilde{\chi}_i &= \chi + \epsilon_{\chi,i} \nonumber \\
\tilde{\eta}_i &= \eta + \epsilon_{\eta,i} \nonumber \\
\tilde{a}_{\alpha,i} &= a_{\alpha} + \epsilon_{a_{\alpha,i}} \nonumber \\
\tilde{a}_{\beta,i} &= a_{\beta} + \epsilon_{a_{\beta,i}} \nonumber
\end{align}
$$

where $\epsilon_{\theta,i}$ with $\theta = \lambda, \varphi, h, \delta, \chi, \eta, a_{\alpha}, a_{\beta}$ are additive noise terms. The propagation of these noisy measurement variables through the observation function $f$ and the knowledge of a reference gravity field allow the estimation of the induced errors $\epsilon_{g,i}$ and $\epsilon_{ggt,i}$ affecting the gravity vector and the gravity gradient tensor respectively given by:

$$
\begin{align}
\epsilon_{g,i} &= g_P - \tilde{g}_{P,i} \nonumber \\
\epsilon_{ggt,i} &= \nabla g_P \cdot u_{\beta\alpha} - \nabla \tilde{g}_{P,i} \cdot u_{\beta\alpha} \tag{17}
\end{align}
$$

By repeating this process, say N times, the estimation of the error expectancy and its variance can be carried out, by means of the following formulas:

$$E[\epsilon_g] = \frac{1}{N} \sum_{i=1}^{N} \epsilon_g, \qquad (18)$$

$$VAR[\epsilon_g] = E[(\epsilon_g - E[\epsilon_g])^2] \qquad (19)$$

$$E[\epsilon_{ggt}] = \frac{1}{N} \sum_{i=1}^{N} \epsilon_{ggt,i} \qquad (20)$$

$$VAR[\epsilon_{ggt}] = E[(\epsilon_{ggt} - E[\epsilon_{ggt}])^2] \qquad (21)$$

Depending on the characteristics of the selected noise models to represent the faults of the measurement system, the magnitude of the gravity field error can then be greater than the expected accuracy. In order to reduce the level of the noise affecting the gravity, a first simple solution consists in averaging the output signal on a sliding window of width $L$, which is equivalent to performing low-pass filtering. The evolution of the variance of the error according to the window size $L$ is then used to characterize the accuracy of the instrument in terms of spatial resolution.

## 4.2 Generation of reference data

To achieve the previously described numerical simulation, it is still necessary to generate reference data consisting of the gravity field related quantities ($g_{P,\alpha,\beta}^n$ and $\nabla g_{P,\alpha,\beta}^n$), submersible trajectories ($\lambda_P, \varphi_P$ and $h_P$), attitude angles ($\delta, \chi$ and $\eta$) and measured accelerations ($a_{\alpha}^b$ and $a_{\beta}^b$).

The reference gravity field derives from a geological model of oceanic crust. The latter combines a bathymetric survey (Figure 4) and an assumed distribution of mineral blocks. Bathymetry is extended to 5000 m depth and the average density of the undersea terrain assumed to be equal to 2.70 g.cm⁻³. The mineral blocks are added to the bathymetric model by introducing local higher densities equal to 3.85 g.cm⁻³. The thickness of the blocks is set at 200 m from the block top surface represented on Figure 5. At this local gravity field is added a global field generated from an ellipsoidal Earth model whose density is set at 5.52 g.cm⁻³. This ensures that the resulting synthetic gravity field has realistic values and variations.

The reference trajectory of the AUV is generated by analyzing actual navigation data which derive from the test mission carried out by the IFREMER off the Mediterranean coasts in the south

<figure><img src="image_5.png" /><figcaption>Figure 4: Perspective view of the bathymetric surface. Width and length of the area are worth about 6 km.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: View of mineralogic block top surface.</figcaption></figure>

of France in 2012. The analysis of these navigation data permits navigation deterministic models to be determined. These models consist of polynomial functions for representing the submersible movement variations at large spatial scales and periodic functions for modeling the variations at smaller spatial scales. Finally, 12 orthogonal profiles were generated. Each profile is of 3600 m long. At an average navigation speed of 4.3 km.h⁻¹, the time travel needed to cover the whole profile is of 3000 s. Note that only the height, $h_P$, is set at −2200 m in all cases. Figures 6 and 7 superimposed the gravity field variations on the 12 profiles on the considered area.

<figure><img src="image_7.png" /><figcaption>Figure 6: Superposition of 12 profiles with variations of gravity field vector components. Results of Part 5. and 6. will be relative to the profile number 3, which has been highlighted on the first map.</figcaption></figure>

Finally, accelerations $a_{\alpha}^b$ and $a_{\beta}^b$ are deduced from the equations of moving-base gravimetry introduced in part 3. Indeed, specific forces in the b-frame can be derived from equations (5) and (6),
<figure><img src="image_8.png" /><figcaption>Figure 7: Superposition of 12 profiles with variations of gravity gradient tensor components.</figcaption></figure>

thus giving:

$$a_{\alpha}^{b} = C_{e}^{b} \ddot{X}_{P}^{c} + C_{b}^{n} (\Omega_{eb}^{b} \Omega_{eb}^{b} + \dot{\Omega}_{eb}^{b}) L_{\alpha}^{b} - C_{n}^{b} g_{\alpha}^{n} \quad (22)$$

and

$$a_{\beta}^{b} = C_{e}^{b} \ddot{X}_{P}^{c} + C_{b}^{n} (\Omega_{eb}^{b} \Omega_{eb}^{b} + \dot{\Omega}_{eb}^{b}) L_{\beta}^{b} - C_{n}^{b} g_{\beta}^{n} \quad (23)$$

Equations (22) and (23) relate the specific forces in the b-frame to reference navigation data and gravity data.

## 4.3 Noise models

As discussed in part 1., the navigation of the Autonomous Underwater Vehicle is achieved through the combination of measures obtained from an Inertial Navigation System, a Doppler Velocity Log and an Ultra Short Base Line acoustic positioning system. The filtering applied to these different observations delivers only the final position and attitude of the carrier and their respective uncertainties. In other words, we have no access to the underlying process leading to the calculation of these quantities. Concerning the positioning of the carrier (longitude $λ_P$ and latitude $φ_P$), the manufacturer indicates an uncertainty equal to 0.1% of the travelled distance. In our case, this represents an uncertainty of about 36 meters at the end of each profile. One possible modelisation of the underlying stochastic process to this position fault can be performed by double integration of a Gaussian noise process. Indeed, in a first approximation, the position delivered can be seen as the result of the double integration of noisy accelerations measured by the inertial unit installed in the INS. The elevation data (or depth $h_P$), will be considered later as being assigned of an uncertainty equal to a few tens of centimeters.

Regarding the attitude angles, the manufacturer indicates an uncertainty about 0.02 deg secant latitude for the heading angle, $δ$, and uncertainties about 0.01 deg secant latitude for the pitch, $χ$, and roll, $η$, angles. Thus, the errors associated with the depth and attitude data are modeled in a first approximation as a Gaussian noise.

On the other hand, faults affecting accelerometers which constitute GRAVIMOB need to be analyzed more carefully. Indeed, in addition to uncertainty assigned to the values they deliver, the sensors are subject to a bias and a scale factor. Although the calibration step aims to set the values of these parameters, they still tend to vary over long time scales. Among the analyzing tools for the study of stochastic process affecting signals, the Allan Deviation is notable for its ease of implementation and its interpretation. Its theory, calculation algorithm and connection with the Power Spectral Density were established in (IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Interferometric Fiber Optic Gyros, 1998). We show the results of the calculation of the Allan Deviation applied to the QFLEX-QA-3000-020 accelerometers (Figure 8). The calculation is performed with three series of static observations of about 8-hour duration at a sampling frequency of 31.25 Hz. According to the properties of Allan Deviation, it is possible to identify three stochastic processes affecting the accelerometer measurements. The first identified process corresponds to a White Noise (WN)

<figure><img src="image_9.png" /><figcaption>Figure 8: Allan Deviation of three QFLEX-QA-3000-020 accelerometers constituting GRAVIMOB. The figure shows clearly the presence of three distinct stochastic processes: white noise, random walk and bias instability.</figcaption></figure>

on the measured accelerations. It is characterized by a slope of $-\frac{1}{2}$ on the Allan Deviation graph. Its standard deviation equal to $7.10^{-6}$ m.s$^{-1} = 7.10^{-1}$ mGal can be read for an averaging time of 1 second on the graph and confirms the manufacturer value. The second process is characterized by a slope equal to $+\frac{1}{2}$ and corresponds to a Random Walk (RW) on the observations. These first two stochastic processes can be easily reproduced numerically from the values given on the graphs of the Allan Deviation. The last process, which is characterized by a zero slope, corresponds to a Bias Instability. Unlike other processes, it is not easy to imitate numerically.

## 5. RESULTS

This section summarizes the results of three numerical tests using reference data and noise models that have been introduced previously (Table 1). For each simulation, 1000 random samples were generated to estimate the standard deviation of the error on the components of the gravity field ($\sqrt{VAR[ε_g]}$ and $\sqrt{VAR[ε_{ggt}]}$). Figures 9, 10 and 11 show the evolution of their maximum value on the profile number 3 (Figure 6) as a function of the width $L$ of the window used in the basic filtering process discussed in part 4.
<table><caption>Table 1: Noise models associated to the numerical simulations. IRW stands for Integrated Random Walk, RW for Random Walk and WN for White Noise.</caption><thead><tr><th></th><th>Test 1</th><th>Test 2</th><th>Test 3</th></tr></thead><tbody><tr><th>Position</th><td>IRW</td><td></td><td></td></tr><tr><th>ε<sub>λ</sub></th><td>IRW</td><td></td><td></td></tr><tr><th>ε<sub>φ</sub></th><td>IRW</td><td></td><td></td></tr><tr><th>ε<sub>h</sub></th><td>WN</td><td></td><td></td></tr><tr><th>Attitude</th><td></td><td>WN</td><td></td></tr><tr><th>ε<sub>δ</sub></th><td></td><td>WN</td><td></td></tr><tr><th>ε<sub>χ</sub></th><td></td><td>WN</td><td></td></tr><tr><th>ε<sub>η</sub></th><td></td><td>WN</td><td></td></tr><tr><th>Accelerations</th><td></td><td></td><td>WN + RW</td></tr><tr><th>ε<sub>a<sub>α</sub></sub></th><td></td><td></td><td>WN + RW</td></tr><tr><th>ε<sub>a<sub>β</sub></sub></th><td></td><td></td><td></td></tr></tbody></table>

## 5.1 Test 1

The first simulation aims to study the effect of the AUV positioning uncertainty. As previously mentioned, one possible modelisation of the underlying stochastic process affecting the position ($ε_{λP}$ and $ε_{φP}$) can be performed by double integration of a White Noise process, or simple Integration of a Random Walk process, abbreviated IRW. The elevation error ($ε_{hP}$) is modeled as a simple White Noise, abbreviated WN.

<figure><img src="image_10.png" /><figcaption>Figure 9: Effect of the AUV positioning uncertainty. At 1 km resolution, the uncertainty affecting $g^U$ is 2 order of magnitude greater than $g^E$ and $g^N$. Components of the gravity gradient tensor are not affected by the AUV positioning uncertainty.</figcaption></figure>

The uncertainties affecting the gravity vector components decrease appreciably when the integration length increases. The graph related to $g^E$ and $g^N$ indicates that the uncertainties affecting these components fall down to 1 mGal at 1 km resolution. However, the uncertainty affecting $g^U$ is 2 order of magnitude greater at the same resolution. It is consistent with the fact that the uncertainty affecting the coordinate $h_P$ is more likely to perturb the restitution of the vertical component $g^U$. The components of the gravity gradient tensor present low uncertainties. This is consistent with equation (13) in which the positioning term does not appear.

## 5.2 Test 2

This simulation focuses on the effect of the uncertainty of the attitude angles ($ε_δ$, $ε_χ$ and $ε_η$). These faults are modeled using a simple White Noise process, WN.

In this case, we still observe the decrease of the uncertainties when the integration length increases. Unlike the previous test, the uncertainty affecting the component $g^U$ falls down to 1 mGal at only 15 m resolution and uncertainties affecting $g^E$ and $g^N$ are 2 order of magnitude greater at the same resolution. Despite an integration over a length of 1 km, uncertainties affecting the components of the gravity gradient tensor range from 100 E and 1000 E.

<figure><img src="image_11.png" /><figcaption>Figure 10: Effect of the AUV attitude uncertainty. Vertical component $g^U$ is less affected than horizontal components $g^E$ and $g^N$. Despite the low-pass filtering, uncertainties affecting the components of the gravity gradient tensor are still high when compared to the gravity effect of geological sources.</figcaption></figure>

## 5.3 Test 3

Here we look at the impact of faults affecting accelerometers ($ε_{aα}$ and $ε_{aβ}$). Allan Deviation led to the identification of different stochastic processes: White Noise, Random Walk and Bias Instability. Because Biais Instability is not easy to imitate numerically, only White Noise, WN and Random Walk, RW processes have been considered in this test.

<figure><img src="image_12.png" /><figcaption>Figure 11: Effect of the accelerometer measurements uncertainty. The low-pass filtering has no effect on the noise reducing because of the non-stationary nature of the first order random walk process involved in the noise model.</figcaption></figure>

Noises on acceleration measurements have a small impact on the gravity vector components. Indeed, the uncertainty is about 1 mGal without applying the low-pass filter. On the contrary, uncertainties on the gravity gradient tensor components are about 60000 E which prevents a proper restitution of the latter. Note that the low-pass filtering here does not reduce noise as in the two previous cases. This can be explained by the non-stationary nature of the first order random walk process introduced into the noise model, which cannot be accounted for by non-predictive filtering procedure.

# 6. KALMAN FILTERING

To improve the estimation, our filtering method must take into account the statistical features of the parameters and noise affecting the observations. The best estimator for this kind of problem is the Kalman Filter (KF). Unfortunately, its original form can not be applied here because of the nonlinearity of the observation equation that links accelerometer measurements to the components of the gravity field. Its extension, the Extended Kalman Filter (EKF) has already been implemented for the data processing of the LiMo-g system (Saint-Jean, 2008). The method allowed the gravity to be estimated from LiMo-g measurements along profiles at a rate equal to the measurement sampling rate. Tests of the method on semi-synthetic data have proved to be conclusive only in short term, that is for short length profiles. Indeed, uncertainty on gravity estimates significantly increases after only
a few stages of calculation, indicating that the Extended Kalman Filter used is likely to be unstable in the long term. One possible way to cope with the problem would be to modify the time evolution model within the Kalman Filter so as to allow for the spatial variability of Earth's gravity field deduced from geostatistics. Theoretically speaking, Markov's autoregressive moving average models can be used for this task (Verdun et al., 2013).

Another complementary solution to the latter is the implementation of the Unscented form of the Kalman Filter (UKF) which was introduced in (Julier and Uhlmann, 1997). The main difference between the EKF and the UKF stems from the manner in which the state probability distribution is propagated through the nonlinear system. In the EKF, the latter is propagated analytically through a first-order linearization while, the UKF uses a minimal set of chosen sample points (called sigma points) which are propagated through the true nonlinear system. It can be shown that the a posterior mean and covariance are captured accurately to second order. On the other hand, the EKF only achieves first-order accuracy. Moreover, the UKF algorithm does not require the determination, implementation and calculation of Jacobian matrix and its complexity is the same order as that of the EKF.

The UKF algorithm is described herein. We use the same notation as in (Haykin, 2001). Let us consider the following discrete-time nonlinear dynamical system:

$$ \begin{cases} x_{k+1} = F(x_k, u_k, v_k) \\ y_k = H(x_k, n_k) \end{cases}, \quad (24) $$

where $x_k$ is the state of the system, in our case the related gravity field quantities that is gravity vector components and gravity gradients, $u_k$ is a known input, $y_k$ is the observed measures of the state, accelerations in our case, $v_k$ is the process noise and $n_k$ is the observation noise.

Like its Extended form, the UKF starts with an initialization phase:

$$ \begin{aligned} \hat{x}_0 &= \mathbb{E}[x_0] \\ P_0 &= \mathbb{E}[(x_0 - \hat{x}_0)(x_0 - \hat{x}_0)^T] \end{aligned} \quad (25) $$

The state vector is then augmented with the process and observation noises:

$$ \begin{aligned} \hat{x}_0^a &= \mathbb{E}[x_0^a] = \begin{bmatrix} x_0^T & 0 & 0 \end{bmatrix}^T \\ P_0^a &= \mathbb{E}[(x_0^a - \hat{x}_0^a)(x_0^a - \hat{x}_0^a)^T] = \begin{bmatrix} P_0 & 0 & 0 \\ 0 & R^v & 0 \\ 0 & 0 & R^n \end{bmatrix} \end{aligned} \quad (26) $$

where $x^a = [x^T v^T n^T]$, $R^v$ is the process-noise covariance and $R^n$ is the measurement-noise covariance. Let $N$ be the length of the augmented state vector $\hat{x}_k^a$. Then, for $k \ge 1, 2N+1$ sigma-points are calculated according to the following deterministic sampling:

$$ \begin{aligned} \chi_{0,k-1}^a &= \hat{x}_{k-1}^a, \\ \chi_{i,k-1}^a &= \hat{x}_{k-1}^a + \left(\sqrt{(N+\lambda)P_{k-1}^a}\right)_i, \quad i = 1, \dots, N \\ \chi_{i,k-1}^a &= \hat{x}_{k-1}^a - \left(\sqrt{(N+\lambda)P_{k-1}^a}\right)_{i-N}, \quad i = N+1, \dots, 2N \end{aligned} \quad (27) $$

where $\lambda$ is a scaling parameter. They can be written as:

$$ \chi_{k-1}^a = \left[ \hat{x}_{k-1}^a \quad \hat{x}_{k-1}^a + \gamma \sqrt{P_{k-1}^a} \quad \hat{x}_{k-1}^a - \gamma \sqrt{P_{k-1}^a} \right] \quad (28) $$

where $\gamma = \sqrt{N+\lambda}$. A more general writing leads to:

$$ \chi^a = \left[ (\chi^x)^T \quad (\chi^v)^T \quad (\chi^n)^T \right]^T \quad (29) $$

These sigma-points are propagated through the nonlinear state function (F):

$$ \chi_{i,k|k-1}^x = F(\chi_{i,k-1}^x, u_{k-1}, \chi_{i,k-1}^v), \quad i = 0, \dots, 2N \quad (30) $$

A Priori estimate and covariance of the state can then be performed as follows:

$$ \hat{x}_k^- = \sum_{i=0}^{2N} W_i^m \chi_{i,k|k-1}^x \quad (31) $$

$$ P_k^- = \sum_{i=0}^{2N} W_i^c (\chi_{i,k|k-1}^x - \hat{x}_k^-) (\chi_{i,k|k-1}^x - \hat{x}_k^-)^T \quad (32) $$

where $W_i^c$ and $W_i^m$ are weights associated to the sigma-points given by:

$$ \begin{aligned} W_0^m &= \frac{\lambda}{N+\lambda} \\ W_0^c &= \frac{\lambda}{N+\lambda} + 1 + \alpha^2 + \beta \\ W_i^m &= W_i^c = \frac{1}{2(N+\lambda)}, \quad i = 1, \dots, 2N \end{aligned} \quad (33) $$

where the constant $\alpha$ controls the spread of the sigma-points and $\beta$ is used to incorporate prior knowledge of the state vector probability distribution. For a Gaussian distribution, $\beta = 2$ is optimal.

A Priori estimation and covariance of the observation are also performed as follows:

$$ \hat{y}_k^- = \sum_{i=0}^{2N} W_i^m \Upsilon_{i,k|k-1} \quad (34) $$

$$ P_{y_k y_k} = \sum_{i=0}^{2N} W_i^c (\Upsilon_{i,k|k-1} - \hat{y}_k^-) (\Upsilon_{i,k|k-1} - \hat{y}_k^-)^T \quad (35) $$

where $\Upsilon_{i,k|k-1}$ are given by:

$$ \Upsilon_{i,k|k-1} = H(\chi_{i,k|k-1}^x, \chi_{k-1}^n), \quad i = 0, \dots, 2N \quad (36) $$

Finally, the Kalman gain matrix is given by:

$$ K_k = P_{x_k y_k} P_{y_k y_k}^{-1} \quad (37) $$

where the cross-correlation matrix $P_{x_k y_k}$ is given by:

$$ P_{x_k y_k} = \sum_{i=0}^{2N} W_i^c (\chi_{i,k|k-1}^x - \hat{x}_k^-) (\Upsilon_{i,k|k-1} - \hat{y}_k^-)^T \quad (38) $$

A Posterior estimate and covariance of the state vector are given by:

$$ \hat{x}_k = \hat{x}_k^- + K_k(y_k - \hat{y}_k^-) \quad (39) $$

and

$$ P_k = P_k^- - K_k P_{y_k y_k} K_k^T \quad (40) $$

We propose here a simple implementation of the Unscented Kalman Filter based on the same synthetic data described in part 4. In this example, only accelerometer observations are assigned a Gaussian noise whose variance is equal to 1 mGal and we seek to restore only the three components of the gravity field vector $g_P^b$. The discrete-time system can be written as:

$$ \left\{ \begin{aligned} g_{P,k+1}^n &= g_{P,k}^n + v_k \\ C_b^n \left( C_e^b \ddot{X}_P^e - a_{P,k}^b \right) &= g_{P,k}^n + n_k \end{aligned} \right. \quad (41) $$

where $a_{P,k}^b = \frac{a_{\alpha,k}^b + a_{\beta,k}^b}{2}$. For $k = 0$, the augmented state vector
and its covariance matrix are given by:

$$ \hat{x}_0^a = [ g_{P,0}^n \quad E[v_0] \quad E[n_0] ]^T \quad (42) $$

$$ P_0^a = \begin{bmatrix} P_0 & 0 & 0 \\ 0 & R^v & 0 \\ 0 & 0 & R^n \end{bmatrix} \quad (43) $$

where $P_0 = \sigma_P^2 I_3$ with $\sigma_P = 0.01 \text{ mGal}$, $R^v = \sigma_v^2 I_3$ with $\sigma_v = 0.1 \text{ mGal}$ and $R^n = \sigma_n^2 I_3$ with $\sigma_n = \frac{1}{\sqrt{2}} \text{ mGal}$.

Figure 12 shows for each component of the gravity field vector, the superposition of the estimate and the reference, the noise level added to the acceleration $a_P^b = \frac{a_{\alpha,k}^b + a_{\beta,k}^b}{2}$ and the noise level on the estimated component which is basically the difference between the true and the estimate value. As previously mentioned, noise level on $a_P^b$ is $\sigma_n = \frac{1}{\sqrt{2}} \text{ mGal}$. The noise level obtained on the gravity field is only about 0.15 mGal which represents a noise level reduction of about 4.5 and demonstrates the efficiency of the basic filter implemented.

<figure><img src="image_13.png" /><figcaption>Figure 12: On a simple example, the Unscented Kalman Filter induces a sensible noise level reduction of about 4.5 on the three components of the gravity vector.</figcaption></figure>

## 7. CONCLUSIONS

We have developed an instrument for measuring the gravity field components in underwater environment. Based on synthetic data, we studied their restitution by introducing different noise processes. The results show that the recovery of the vertical gravity vector component with an accuracy suitable to geological sources detection is feasible by improving the positioning of the AUV. The recovery of gravity vector horizontal components and the gravity gradients requires very good estimates of carrier attitude angles, and thus high grade attitude and heading reference sys-

tems are needed. It must be emphasized that whatever the performances of positioning and attitude determination systems, the uncertainty affecting the gravity related quantities depends directly on the accelerometer sensibility. Given one type of accelerometers, the unique way to reduce uncertainties on gravity estimates relies upon the use of predictive filtering methods able to account for non-stationary noise such as random walk which affects this kind of sensor. In this work, an Unscented Kalman Filter has been implemented in a simple case. The continuation of this work will be the inclusion of all stochastic processes affecting the measures and the covariance of the gravity field components. These new refinements will make possible even more accurate undersea gravity measurements.

## ACKNOWLEDGEMENTS

The authors are indebted to the French Ministry of Defence and the Pays de la Loire Region for their support of this work.

## REFERENCES

Araya, A., Kanazawa, T., Shinohara, M., Yamada, T., Fujimoto, H., Iizasa, K. and Ishihara, T., 2011. A gravity gradiometer to search for submarine ore deposits. In: Underwater Technology (UT), 2011 IEEE Symposium on and 2011 Workshop on Scientific Use of Submarine Cables and Related Technologies (SSC), Tokyo, Japan, pp. 1–3.

Haykin, S. S., 2001. Kalman filtering and neural networks. Adaptive and learning systems for signal processing, communications, and control, Wiley, New York.

IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Interferometric Fiber Optic Gyros, 1998. IEEE Std 952-1997.

Julier, S. J. and Uhlmann, J. K., 1997. New extension of the Kalman filter to nonlinear systems. In: Signal Processing, Sensor Fusion, and Target Recognition VI, Vol. 3068, Ivan Kadar, Orlando, Florida, USA, pp. 182–193.

Liu, F., Qian, D., Zhang, Y. and Li, Y., 2010. A computer simulation of the influence of GGI and inertial sensors on gravity gradient aided navigation. In: 2010 3rd International Symposium on Systems and Control in Aeronautics and Astronautics (ISSCAA), Harbin, China, pp. 793–797.

Saint-Jean, B. D., 2008. Étude et développement d'un système de gravimétrie mobile. PhD thesis, Observatoire de Paris.

Verdun, J., Damenet, N. and Cali, J., 2013. Moving-base vector gravimetry data processing based on optimal, physically sensible evolution models. In: Proceedings of International Symposium, Terrestrial Gravimetry: Static and Mobile Measurements, Elektropribor International Association of Geodesy, TG-SMM 2013, St Petersburg, Russia, pp. 41–47.

Yan, Z., Yang, F., Ma, J. and Tian, J., 2012. Underwater Obstacle Detection Based on the Change of Gravity Gradient. In: 2012 2nd International Conference on Remote Sensing, Environment and Transportation Engineering (RSETE), Jiangsu, China, pp. 1–5.

Yao, Z. and Xiaorong, S., 2012. Densification algorithm research on gravity gradiometer reference map generation. In: 2012 International Conference on System Science and Engineering (IC-SSE), Coventry, United Kingdom, pp. 315–320.

Revised March 2015