Next Article in Journal
An Investigation on Super- and Sub-Terminal Drops in Two Different Rain Categories and Climate Regimes
Previous Article in Journal
Evaluation of Sea Ice Radiative Forcing according to Surface Albedo and Skin Temperature over the Arctic from 1982–2015
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development of a Lightweight Inertial Gravimeter for Use on Board an Autonomous Underwater Vehicle: Measurement Principle, System Design and Sea Trial Mission

1
Laboratoire Géomatique et Foncier (GeF) EA 4630, Conservatoire National des Arts et Métiers (Cnam), HESAM Université, École Supérieure d’Ingénieurs Géomètres Topographes (ESGT), 1 Boulevard Pythagore, 72000 Le Mans, France
2
Laboratoire Géosciences Océan (LGO)—UMR6538, Institut Universitaire Européen de la Mer, IUEM—Technopôle Brest Iroise, rue Dumont D’Urville, 29280 Plouzané, France
3
Mappem Geophysics, Bâtiment Tech-Iroise, 1 rue des Ateliers, Zone de Mespaol, 29290 Saint-Renan, France
4
Service Positionnement, Robotique, Acoustique et Optique (PDG-DFO-SM-PRAO), IFREMER Centre Méditerranée, Zone Portuaire de Brégaillon, CEDEX CS20 330, 83507 La Seyne-sur-Mer, France
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(11), 2513; https://doi.org/10.3390/rs14112513
Submission received: 31 March 2022 / Revised: 15 May 2022 / Accepted: 16 May 2022 / Published: 24 May 2022
(This article belongs to the Section Engineering Remote Sensing)

Abstract

:
The purpose of this paper is to present the design, development and testing of an innovative instrument called GraviMob, which allows performing dynamic measurements of underwater gravity anomalies. After recalling the interest in underwater gravimetry, we describe the system, the core of which consists of triads of accelerometers rigidly attached to an Autonomous Underwater Vehicle (AUV). The article also presents the mathematical methods for estimating the east, north and vertical components of the local gravity vector. An unscented Kalman filter, integrating AUV position and orientation data, performs estimation of gravity in a frame adapted to its interpretation. To assess its performance, GraviMob was tested in the Mediterranean Sea during the year 2016. A comparison of the surface gravimetric signal previously acquired by the French Navy indicates that the maximum discrepancy between the vertical gravity component and its reference is below 4 mGal. Components of the vertical deflection calculated from GraviMob’s measurements were compared with those calculated from recent gravity field models. While a remarkable agreement was found on the north component, there remains a discrepancy (7 arcsec) on the east component which can be largely reduced by refining the estimation of the orientation of GraviMob’s sensitive axes in the AUV.

Graphical Abstract

1. Introduction

Deep-sea exploration and sustainable resource exploitation have generated considerable research interest in both basic and applied geophysical research. The development of new technologies for deep-sea vehicles such as high-resolution acoustic sonars, multibeam echo-sounders and magnetometers, allows very detailed images of the seafloor to be produced. These developments have not only improved our understanding of the dynamics of major faults, mid-oceanic ridges and hydrothermal vents but also opened new perspectives on the exploration and exploitation of mineral resources in deep water [1,2,3]. Today, the study of the structure of the seafloor at tens to a few kilometres largely relies on high-resolution bathymetry and magnetism, particularly when conducted close to the seafloor onboard Autonomous Underwater Vehicles (AUVs). Meanwhile, sea surface gravity surveys carried out onboard oceanographic vessels have revolutionised our vision of the structure of the oceanic crust revealing a high spatial and temporal variability [4,5,6]. However, the water depth that separates gravity measurements from gravity sources (2 to 5 km) greatly attenuates most high-frequency gravity signals, indispensable to imaging the smallest seafloor geological structures. Some rather old experiments of seafloor deployment of a customised static gravimeter installable on sea-bottom from a submarine or a Remotely Operated Vehicle (ROV), have significantly increased our knowledge of the paths of percolation of seawater and fracturing [7] and of the in-depth distribution and density of hydrothermal sulphide deposits [8]. Despite these indisputable advances, such experiments have proven to be less efficient in terms of spatial resolution since the benefit of manoeuvring near seafloor gravity sources is offset by the uneven spatial distribution of gravity measurements. Today, available mobile gravimetric systems are still cumbersome and energy-intensive, which in particular prohibits their installation on light submarine vehicles, such as AUVs. Hence, the development of a new type of gravimetric sensor with a small footprint and low energy consumption, suitable for installation onboard AUV appears essential to meet the problem posed by the measurement of the Earth’s gravity field near seafloor with unprecedented accuracy and spatial resolution.
The development of gravity sensors for underwater vehicles is a key project for some of the most prestigious institutes working on ocean sciences, most notably the Woods Hole Oceanographic Institute, USA, and the Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Japan. Research at JAMSTEC and the University of Tokyo is probably the most advanced. The Japanese system consists of a LaCoste-Romberg S-174 Micro-g gravimeter mounted on a gimbal mechanism and a gradiometer consisting of two astatic pendulums also mounted on a gimbal platform [9,10,11]. Since 2012, this system has been repeatedly tested on board the submersible URASHIMA [12,13,14,15]. On the basis of sea trials, the accuracy of this system is of the order of a tenth of mGal (1 mGal = 10 5 m s 2 ) for gravimetric measurement (RMS of cross over errors after levelling) [16]. To our knowledge, the Japanese team is the first to have calculated and interpreted Bouguer gravity anomalies from AUV data acquired at sea. The dimensions of the Japanese AUV (1.3 m in width, 10 m in length, 1.5 m in height), much larger than those of the other AUVs, allowed them to adopt a classical approach, with a sensor based on a terrestrial gravimeter (a Lacoste-Romberg) mounted on a compensation system controlled by gyros. Their prototype is, however, very greedy in energy and space and difficult to transpose to other AUVs of smaller dimensions.
For its part, the American team from the Department of Geology and Geophysics at the Woods Hole Oceanographic Institute is focusing its efforts on modelling the navigation of autonomous underwater vehicles [17,18,19]. On the basis of numerical simulations, the team is seeking to reduce, by filtering and smoothing, the impact of uncertainties affecting the vertical positioning of the underwater vehicle on the restitution of the underwater gravity field. This work does not mention the development of new sensors, which suggests the adaptation and use of terrestrial gravity instruments, such as relative gravimeters, on board autonomous underwater vehicles that would be verticalised by means of a stabilised platform.
A recent experiment of gravity measurement onboard an unmanned underwater vehicle (UUV) and using a strapdown sensor was carried out by a Chinese team in Mulan Lake, Wuhan City, Hubei Province, China [20]. The survey was carried out using a small UUV BQR800 (1.55 m in length, 0.87 m in width, 0.81 m in height) remotely operated at a depth of 8 to 20 m along a 6 km profile. The gravity measurements were delivered by a dg-M shipborne strapdown gravimeter of 1 mGal accuracy developed by Hunan INS Technology Co., Ltd. Submersible positioning and orientation were ensured by a set of positioning sensors consisting of an inertial navigation system, a short baseline acoustic system, (accuracy: 1 m at 100 m distance) an ultrasonic Doppler velocimeter (accuracy: 0.1 cm/s) and a depth gauge (accuracy: 2 mm). According to the Chinese team’s findings, based on the analysis of six repeated surface and underwater profiles, the accuracy of the system reached 0.42 mGal. These results tend to demonstrate the feasibility of mobile underwater gravimetry only in lakes and need to be confirmed in the underwater environment at greater depths, where navigation conditions are more difficult.
The purpose of this study is to present the design, development and testing of an innovative mobile gravity sensor that allows the dynamic measurement of the Earth’s gravity anomalies near the ocean floor. This instrument, named GraviMob (Mobile Gravimetry System), is inspired by the LIMO-g system (Light Mobile Gravimetry System) whose development was initiated in 1999 by the French National Institute of Geographic and Forest Information (IGN), then studied successively during the doctorates of Bertrand de Saint-Jean [21] and Clément Roussel [22]. Compared to competing systems, GraviMob is the only one that can measure the gravity vector by determining its east, north and vertical components. GraviMob is a lightweight, low power consumption sensor that can be easily integrated on a small AUV such as those used by the French Oceanographic Fleet (FOF). In this article, Section 2 is entirely devoted to the description of the GraviMob system and the data processing methods. We first present the sensor we have developed, and the carrier vehicle used to test the GraviMob system (Section 2.1). The mathematical model of the GraviMob system (Section 2.2) and the methodology used to determine the internal parameters of the system (Section 2.3) are then discussed in detail. Finally, the method of processing the GraviMob measurements together with the positioning and orientation data of the carrier vehicle is fully explained (Section 2.4). The results of the GraviMob system experimentation in the Mediterranean Sea are described in Section 3. Planned modifications to the GraviMob system to improve its performance are discussed in Section 4. Additional information necessary for the understanding of the article is given in the appendices. Definition of useful coordinate reference frames and formulae of coordinate changes are given in Appendices Appendix A and Appendix B. Explicit forms of transformation and skew symmetric matrices are provided in Appendix C. Finally, a constrained optimisation method is described in Appendix D.

2. Materials and Methods

2.1. GraviMob System Composition

GraviMob system essentially consists of one accelerometer sensor which itself contains six QFLEX-QA 3000 unidirectional electrostatic accelerometers manufactured by the American company Honeywell (Figure 1a). The working principle of each electrostatic accelerometer is the following: a proof mass hanging to a flexible beam and exposed to a magnetic field from an electromagnet, is kept at rest thanks to the restoring magnetic force exerted on it. If the acceleration of the proof mass varies, the resulting displacement of the proof mass relative to the sensor box is measured by a capacitive displacement sensor, that generates a signal proportional to the displacement amplitude. This signal is then fed back to the control system that determines the intensity of the current needed to supply the electromagnet so as to make the proof mass return to its rest position.
The six accelerometers are mounted in two triads whose two-by-two orthogonal axes allow the components of the specific force vector to be measured (Figure 1b). Moreover, these two triads are rigidly mounted in the carrier vehicle used for performing gravity surveys (Figure 1c,d). Each accelerometer from one of the two triads provides a measurement of the specific force component along its sensitive axis. The specific force is tantamount to the force per unit mass that has to be exerted on the proof mass of one given accelerometer to maintain it motionless with respect to the sensor box. By itself, the sole measurements of the specific force components do not permit those of the gravitational force to be estimated. They actually need to be corrected from the motion-induced accelerations experienced by the proof mass. For that reason, GraviMob jointly makes use of position data provided by an additional positioning system. Such position data are commonly expressed in a terrestrial reference frame, whose axes are not aligned with accelerometer sensitive axes. This is why, as we shall see subsequently through mathematical relationships (see Section 2.2), the joint determination of the carrier vehicle orientation is essential. In our experiments, the Autonomous Underwater Vehicle (AUV) Aster x , operated by IFREMER, was used as the carrier vehicle (Figure 2).
A sophisticated navigation system on board allows measurements of the position and the orientation of the AUV to be continuously acquired. The navigation system consists of three distinct devices: an inertial navigation system (INS) from iXBlue company, an ultrasonic Doppler velocimeter and an ultra-short baseline (USBL) positioning system, that calculates the position by measuring the range from a vessel-mounted transceiver to an acoustic transponder fitted to the AUV. The vessel in question, which also supports the AUV operation by putting and getting back the submersible to its working site, is itself accurately positioned with the aid of GNSS techniques. Unlike shipborne or airborne gravity measurements, those acquired in the submarine domain are relatively undisturbed by submersible-induced accelerations. Consequently, the measuring range of the accelerometers can be significantly narrowed (Table 1), thus lessening the digitizing error at 0.02 mGal using a 24 bit digitiser with a sampling frequency at 2 Hz .
The accelerometers’ noise level proved to be at 1 mGal / Hz based on an experimental noise analysis drawn upon Allan variance [23]. Through the combined processing of all positioning sensors’ data, the overall planimetric and altimetric accuracies of AUV positioning have been estimated at 1 to 2.5 m and 30 cm respectively. The attitude angles jointly determined consist of the yaw angle measured at 0.05 deg accuracy and the pitch and roll angles both measured at 0.005 deg accuracy. The whole apparatus including the two accelerometer triads and the 24 bit digitiser, is enclosed in a 17 diameter ( 43.2 cm ) watertight sealed glass sphere. The sphere is locked in an ad hoc plastic frame that is screw-fixed into the front case of the AUV. The overall weight of the system is of some 40 kg and its proper function requires a power supply of 2.5 W. It should be noted that the sphere is entirely immersed during gravity survey operations, thus allowing a thermal balance between the gas inside the sphere and sea water to be achieved. The temperature stability inside the sphere stems exclusively from this thermal balance, with no need for an additional temperature-stabilizing system.

2.2. GraviMob System Mathematical Model

This subsection is devoted to deriving the fundamental equation that links gravity to the measurements performed by the GraviMob system. We first need to define the coordinate reference frames used in the paper (see Appendix A) and describe some properties of the transformation matrix between two reference frames (see Appendix B). The explicit forms of both transformation matrices and skew symmetric matrices are also given in Appendix C. We draw heavily on the relationships in Jekeli [24], Rogers [25] and Titterton et al. [26] where greater mathematical detail is given. In the following, we will denote P the point of the carrier vehicle at which gravity is actually determined.
Let us consider one 3D accelerometer labelled by α . Let X M α i be the one column matrix which contains the coordinate in the i-frame of the position vector of the proof mass inside the accelerometer case. The proof mass is located at the point M α where acceleration measurement is performed, which also coincides with the origin of the s-frame s α defined by the accelerometer α . According to Newton’s Second Law applied to the proof mass, the second-order time derivative X ¨ M α i of X M α i may be expressed as
X ¨ M α i = g M α i + a M α i ,
where g M α i is the gravitational acceleration and a M α 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 vector X M α i can be related to the position vector of the point P expressed in the e-frame, X P e , as
X M α i = C e i X P e + C b i L α b ,
where C e i and C b i are the direction cosine matrices of the i-frame relative to the e-frame and b-frame respectively, and L α b the lever arm between the point P and the measuring point M α of the accelerometer α . It must be emphasised that the restoring force a M α 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 α b are constants independent of time.
The specific force a M α i is directly related to the accelerometer measurements, that is a M α s α components where s α is the s-frame associated with the accelerometer α , by
a M α i = C n i C b n C s α b a M α s α ,
where C n i and C b n are the orthogonal transformation matrices from the i-frame to the navigation frame n-frame, and from the body frame b-frame to the n-frame respectively, and C s α b is the transformation matrix not necessarily orthogonal from the s-frame to the b-frame. The elements of matrix C s α b can be determined once and for all for each 3D accelerometer during the design of the system. Some authors [27] suggest that these elements could vary in the long term and should therefore be steadily upgraded by means of pre-calibration or on-the-fly methods. The description of such a method will be dealt with specifically in Section 2.3.
By differentiating Equation (2) twice with respect to time according to Equation (A12) and given the fact that L α b does not depend on time ( L ˙ M α b = L ¨ M α b = 0 ), we obtain
X ¨ M α i = C e i X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e + Ω ˙ i e e X P e + C b i Ω i b b Ω i b b + Ω ˙ i b b L α b ,
where Ω i e e is the skew symmetric matrix associated with the rotation of the e-frame with respect to the i-frame and Ω ˙ i e e its first-order time derivative, both expressed in the e-frame, and similarly Ω i b b is the skew symmetric matrix associated with the rotation of the b-frame with respect to the i-frame and Ω ˙ i b b its first-order time derivative, both expressed in the b-frame.
According to Equations (1) and (4), the gravitational acceleration g M α i given in the i-frame may be expressed as
g M α i = C e i X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e + Ω ˙ i e e X P e + C b i Ω i b b Ω i b b + Ω ˙ i b b L α b a M α i .
By multiplying both sides of Equation (5) by matrix C i n using equations C i n C e i = C e n , C i n C b i = C b n , C i n C n i = I 3 and Equation (3), we obtain
g M α n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e + Ω ˙ i e e X P e + C b n Ω i b b Ω i b b + Ω ˙ i b b L α b C b n C s α b a M α s α .
This equation relates the components of the gravitational acceleration to the three vectors corresponding respectively to the position, velocity and acceleration of the point P, the specific force and the lever-arm, all of those being defined in the b-frame. By assuming L α b = 0 in Equation (6) and expressing X ¨ P e as a function of the remaining terms in (6), we obtain the so-called navigation equation in the e-frame already derived by Jekeli [24].
The gravity vector denoted by g ¯ M α n is then deduced from the gravitation vector by subtracting the centripetal acceleration Ω i e e Ω i e e X M α e due to the Earth’s rotation expressed in the n-frame, thus giving:
g ¯ M α n = g M α n C e n Ω i e e Ω i e e X M α e ,
where the direction cosine matrix C e n depends only on the position of point P relative to the e-frame, Ω i e e is the skew symmetric matrix of the Earth’s own rotation around the polar axis and the position vector of point P in the e-frame, X M α e , may be expressed as:
X M α e = X P e + C n e C b n L α b .
The equation giving the gravity vector can be deduced by combining Equations (6)–(8); it thus follows:
g ¯ M α n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω ˙ i e e X P e + C b n Ω i b b Ω i b b + Ω ˙ i b b C e n Ω i e e Ω i e e C n e C b n L α b C b n C s α b a M α s α .
Equation (9) may also be expressed for the second accelerometer triad indicated by β as follows:
g ¯ M β n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω ˙ i e e X P e + C b n Ω i b b Ω i b b + Ω ˙ i b b C e n Ω i e e Ω i e e C n e C b n L β b C b n C s β b a M β s β .
The presence of two accelerometric triads makes it possible to avoid the term relative to the lever arm. This remains true as long as the point P above mentioned is such that P M α = L α = P M β = L β , i.e., in the case where point P is exactly located in the middle of the segment [ M α M β ] . In that case, combining Equations (9) and (10) as an arithmetic average leads to a simple equation where all terms related to the lever arm have cancelled each other out:
1 2 g ¯ M α n + g ¯ M β n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω ˙ i e e X P e 1 2 C b n C s α b a M α s α + C s β b a M β s β .
Under the above hypotheses, Equation (11) can also be regarded as a second-order approximate of the gravity vector g ¯ P n at point P expressed in the n-frame. Indeed, we can express g ¯ M α n and g ¯ M β n in the following expansions about, respectively, the lever arms L α n and L β n :
g ¯ M α n = g ¯ P n + J P X P e L α n + L α n ε α L α n ,
g ¯ M β n = g ¯ P n + J P X P e L β n + L β n ε β L β n ,
where J P X P e corresponds to the Jacobian matrix at point P of the function defined from R 3 to R 3 that relates every position vector X M e to the gravity vector g ¯ M n at point M, for a = α , β , L a n denotes the Euclidian norm of vector L a n , and ε a stands for a vectorial function such that:
ε a L a n L a n 0 0 .
Provided that L α n = L β n , the arithmetic average (11) obtained by Equations (12) and (13) is then given by:
1 2 g ¯ M α n + g ¯ M β n = g ¯ P n + 1 2 L α n ε α L α n + ε β L β n ,
which indicates that the gravity vector g ¯ P n differs from the arithmetic average (Equation (11)) solely by one second order term proportional to L a n 2 . This approximation is all the more accurate as the length of the lever arm L a n is small compared to that of the position vector X P e . We will now adopt the following expression of the gravity vector at point P:
g ¯ P n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω ˙ i e e X P e 1 2 C b n C s α b a M α s α + C s β b a M β s β .
By a similar approach, an expression of gravitational acceleration g P n at point P can be established and expressed as follows:
g P n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e + Ω ˙ i e e X P e 1 2 C b n C s α b a M α s α + C s β b a M β s β .
Among the terms that compose Equation (16) (or (17)) is Ω ˙ i e e X P e , which explicitly depends on fluctuations in the Earth’s angular rate of rotation. Let ω e be the value of the Earth’s angular rate of rotation, then we can state that:
Ω ˙ i e e X P e | ω ˙ e | X P e .
To estimate an upper bound of this term, we consider a major contribution to the time variations of the Earth’s angular rate of rotation, that comes from the bimonthly fluctuation of the Length Of Day ( L O D ) estimated at 0.5 milliseconds. If δ L O D denotes the variation of the length of day during the time interval δ t , then the time variation ω ˙ e of the Earth’s angular rate of rotation may be expressed as:
ω ˙ e = ω e L O D δ L O D δ t .
From the following numerical values [28,29]:
ω e = 7.292115 × 10 5 rad s 1 , L O D = 86 , 400 s , δ L O D = 1 ms , δ t = 0.5 month , X P e = 6 , 378 , 137 m ,
we can obtain the following reasonable estimates of | ω ˙ e | and | ω ˙ e | X P e :
| ω ˙ e | 10 18 rad s 2 ; | ω ˙ e | X P e 10 11 m s 2 = 10 6 mGal .
According to the inequality (18) and given the level of noise of GraviMob accelerometers at 1 mGal / Hz , the term Ω ˙ i e e X P e can be neglected in Equations (16) and (17). Once the term Ω ˙ i e e X P e has been removed, Equation (16) (resp. (17)) is the desired gravity (resp. gravitational acceleration) observation equation:
g ¯ P n = C e n X ¨ P e + 2 Ω i e e X ˙ P e 1 2 C b n C s α b a M α s α + C s β b a M β s β ,
resp . g P n = C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e X P e 1 2 C b n C s α b a M α s α + C s β b a M β s β .
Finally, it should be noted that since both accelerometer triads are placed at the front of the AUV, there is an offset L C b = P C b , where P C b stands for the components of vector P C in the b-frame, between point P and the point C (AUV reference point) effectively positioned by the positioning system of the AUV. The position of point P, X P e , can only be determined providing the position X C e given by the AUV positioning system and the vector L C b joining point P to the point C. The components of this vector can be accurately determined from a 3D geometrical model describing the AUV and the accelerometer triads. The position vector X P e can then be related to the position vector of the point C, X C e , as
X P e = X C e + C n e C b n L C b .
A major improvement in the design of the GraviMob system would be to place the positioning system so that the points P and C coincide ( L C b = 0 ). This configuration would be possible by placing the INS in a central position between the two triads, which could be considered for a future version of the GraviMob system.

2.3. GraviMob System Adjustment

2.3.1. Basic Equation of Adjustment

In practice, the signals delivered by each triad of accelerometers consist of three electrical voltages, which are also affected by biases, i.e., offset voltages even when the sensors are not accelerated. Let V s α be the column vector whose components correspond to the three electrical voltages measured in the triad α . The following reasoning will remain valid in all respects for the triad β . Assuming a linear response of the accelerometers, the relationship that relates this vector to the specific force a M α s α may be expressed as:
V s α = K α a M α s α + V 0 s α ,
where K α is a 3 × 3 diagonal matrix consisting of scale factors and V 0 s α is the vector of the three electrical voltage biases. For the purpose of performing adjustment, the relation (23) has to be slightly modified as follows:
V s α = K α C s o s α a M α s o + V 0 s α ,
where C s o s α is the transformation matrix from the s o -frame to the s-frame corresponding to the triad α . It should be recalled that the axes of s-frame are not necessarily orthogonal due to misalignments of accelerometer sensitive axes. The matrix C s o s α is thus non orthogonal. Moreover, the axes of s o -frame are not necessarily aligned with those of b-frame. This point will be dealt with at the end of Section 2.3. Adjustment of GraviMob system consists of each triad in estimating GraviMob’s internal parameters, that is the bias vector V 0 s α and the elements of matrices K α and C s o s α through experimental painstaking techniques, that can be carried out in laboratory just before the survey.

2.3.2. Overview of an User-Friendly Adjustment Method

In general terms the idea of adjustment methods consists in collecting several measurements of V s α and a M α s α while keeping the system stationary. In that case, the system measures only the gravity vector and thus the amplitudes of the three signals delivered by the triad depend essentially on the gravity intensity and the orientation of the system. Mathematically speaking, if g ¯ M α n denotes the gravity vector at the location of the triad, the specific force a M α s α may be related to the gravity vector given the transformation matrix C n s α from the n-frame to the s-frame by:
a M α s = C n s g ¯ M α n .
By carrying out static measurements with the system tilted in N different orientations (Figure 3), N pairs
a M α s α , V s α i , i = 1 , 2 , , N
can be formed, that are theoretically linked by Equation (24). The system adjustment, therefore, amounts to estimating the best parameter values consistent with the N pairs obtained from static measurements.
Let k x , k y , k z be the three scale factors such that the matrix K α may be expressed as:
K α = k x 0 0 0 k y 0 0 0 k z .
In the case where k x = k y = k z and the reference frames s-frame and s o -frame are orthogonal, the system adjustment is tantamount to solving the so-called orthogonal Procrustes problem whose solution can be found in Schönemann [30]. In the more general case where k x k y k z and s-frame and s o -frame are not orthogonal, the Procrustes solution is no longer valid. It is then necessary to use an iterative search of GraviMob’s internal parameters [31]. The method we adopted is detailed in Yang et al. [32]. It assumes only the prior determination of gravity vector magnitude at the site of experimentation. In order to give the explicit mathematical formulation of the method, the following explicit notations should be used:
V s α k = v x , i v y , i v z , i ,
where v x , i , v y , i , v z , i are the three electrical voltages given by the triad α for the i-th orientation of GraviMob system;
a M α s o k = g ¯ M α n c x , i c y , i c z , i ,
where g ¯ M α n denotes the magnitude of the local gravity vector and c x , i , c y , i , c z , i are the three direction cosines such that g ¯ M α n c x , i , g ¯ M α n c y , i , g ¯ M α n c z , i correspond respectively to the three components of the gravity vector g ¯ M α n once projected onto x, y and z-axes of the s o -frame for the i-th orientation of GraviMob system; this equation is equivalent with Equation (25);
C s o s α = 1 τ x y τ x z 0 τ y y τ y z 0 0 τ z z ,
denotes the non-orthogonal transformation matrix from the s o -frame to the s-frame defined in Appendix C from a formula given in Panahandeh et al. [33];
V 0 s α = v x 0 v y 0 v z 0 ,
where v x 0 , v y 0 , v z 0 are the three electrical voltage biases of the triad α . For the i-th orientation of GraviMob system, Equation (24) may be explicitly expressed as:
v x , i v y , i v z , i = g ¯ M α n k x 0 0 0 k y 0 0 0 k z 1 τ x y τ x z 0 τ y y τ y z 0 0 τ z z c x , i c y , i c z , i + v x 0 v y 0 v z 0 = g ¯ M α n k x x k x y k x z 0 k y y k y z 0 0 k z z c x , i c y , i c z , i + v x 0 v y 0 v z 0 ,
after having defined k l m , l , m = x , y , z terms as follows:
k x x = k x , k x y = k x τ x y , k x z = k x τ x z , k y y = k y τ y y , k y z = k y τ y z , k z z = k z τ z z .
Given an initial set of N distinct orientations ( N > 3 ) defined by direction cosines c x , i , c y , i , c z , i , i = 1 , 2 , , N , the six parameters k x x , k x y , k x z , k y y , k y z and k z z and the three voltage biases v x 0 , v y 0 and v z 0 can be estimated by solving the overdetermined linear system of 3 N equations in 9 unknowns by least squares method expressed as:
A 1 A 2 A 3 A i A N 1 A N k x x k x y k x z v x 0 k y y k y z v y 0 k z z v z 0 = b 1 b 2 b 3 b i b N 1 b N ,
where
A i = g ¯ M α n c x , i g ¯ M α n c y , i g ¯ M α n c z , i 1 0 0 0 0 0 0 0 0 0 g ¯ M α n c y , i g ¯ M α n c z , i 1 0 0 0 0 0 0 0 0 0 g ¯ M α n c z , i 1 ,
and
b i = v x , i v y , i v z , i .
Once those nine parameters have been determined, the 3 N values of the direction cosines c x , i , c y , i and c z , i , i = 1 , 2 , , N can be updated by finding out the minimum of the N cost functions defined as:
f i ( c x , i , c y , i , c z , i ) = K α C s o s α a M α s o + V 0 s α V s α i 2
that can explicitly written as follows
f i ( c x , i , c y , i , c z , i ) = g ¯ M α n k x x c x , i + g ¯ M α n k x y c y , i + g ¯ M α n k x z c z , i + v x 0 v x , i 2 + g ¯ M α n k y y c y , i + g ¯ M α n k y z c z , i + v y 0 v y , i 2 + g ¯ M α n k z z c z , i + v z 0 v z , i 2 .
Assuming that the experiments carried out to tilt the GraviMob system along different directions are independent, the direction cosines c x , i , c y , i and c z , i for each inclination can then be considered statistically independent. They can be thus estimated by minimizing N independent cost functions as (34) under the constraints stating that c x , i 2 + c y , i 2 + c z , i 2 = 1 for all i varying from 1 to N. This minimisation under non-linear constraints is explained in detail in Appendix D. The updated values of the direction cosines are then used to reiterate an estimation of the above mentioned internal parameters and such a recursive process is then resumed until a stopping criterion will be satisfied. To this end, a residual r i ( m ) is calculated for the i-th orientation and the m-th iteration, which may be expressed as:
r i ( m ) = C s α s o K α 1 V s α V 0 s α i g ¯ M α n .
Let σ v ( m ) be the standard deviation of the N residuals r i ( m ) , i = 1 , 2 , , N calculated for the m-iteration. Then, the stopping criterion consists in verifying the following inequality:
| σ v ( m ) σ v ( m 1 ) | < ε ,
where ε is a threshold value given a priori. Such an inequality ensures that the dispersion of the residuals has not changed significantly between the ( m 1 ) -th and the m-th iteration. Once the six parameters defined in Equation (32) have been determined, the three scale factors k x , k y , k z and the six parameters τ x y , τ x z , τ y y , τ y z and τ z z can be calculated provided that three additional relationships involving them can be found. Noting that the columns of the matrix C s o s α 1 = C s α s o are the coordinates of the three non orthogonal unit vectors of s-frame-axes expressed in the s o -frame, we have:
C s α s o = 1 τ x y τ y y τ x y τ y z τ x z τ y y τ y y τ z z 0 1 τ y y τ y z τ y y τ z z 0 0 1 τ z z ,
and thus,
τ y y 2 = 1 + τ x y 2 , τ z z 2 = 1 + τ x z τ y z τ x y τ y y 2 + τ y z τ y y 2 .
Since the misalignment angles between the respective y- and z-axes of the s-frame and s o -frame are infinitesimal, we can infer that the two last diagonal terms of matrix C s α s o are close to 1, and thus positive. As a result, the nine internal parameters of GraviMob can finally be calculated using the following equations used successively:
k x = k x x , τ x y = k x y k x x , τ x z = k x z k x x , τ y y = 1 + τ x y 2 , k y = k y z τ y y , τ y z = k y z k y , τ z z = 1 + τ x z τ y z τ x y τ y y 2 + τ y z τ y y 2 , k z = k z z τ z z .

2.3.3. Testing of the Adjustment Method

After being validated by numerical simulations, the adjustment method was applied to experimental measurements acquired in laboratory, and consisted of the six voltages v k q , k = x , y , z , q = α , β delivered by the accelerometers of the two triads held fixed for different orientations of the GraviMob system. The value of the local gravity intensity in our laboratory is known, equal to 980,856.4908 ± 0.0034 mGal, measured by an FG5 absolute gravimeter (FG5#206) between 26 and 29 April 1999. In order to limit the effect of temperature variation, the acquisition time was deliberately limited to 15 min during which N = 50 different orientations could be tested. The temperature variation observed during the experiment was 0.22 °C for an average temperature of 23.82 °C. The iterative parameter estimation algorithm of the adjustment method was initialised with the initial values of the direction cosines calculated by:
c x , i v x , i v x , i 2 + v y , i 2 + v z , i 2 , c y , i v y , i v x , i 2 + v y , i 2 + v z , i 2 , c z , i v x , i v z , i 2 + v y , i 2 + v z , i 2 ,
for all i varying from 1 to N. Given a threshold value of ε = 10 16 (Equation (36)), the number of iterations necessary for the convergence of the algorithm was approximately 10 , 000 . However, in order to observe the dispersion of internal parameter values, this number was voluntarily fixed at 50 , 000 so as to evaluate the mean and standard deviation of the estimated values of internal parameters. These values obtained with the measurements acquired during this experiment are given in Table 2. The corresponding residuals (Equation (35)) are zero mean. Their standard deviations for the triads α and β are equal to 1.4 and 3.1 mGal respectively. These statistical indicators make it possible to evaluate the uncertainty of GraviMob measurements resulting specifically from the uncertainty of the values of the internal parameters. According to our findings, this uncertainty cannot be less than 3 mGal due to the measurement range of the GraviMob system, for which the possible tilt angles must be strictly within the range of ±10.8 degrees (see Table 1 in Section 2.1). Of particular note is the high uncertainty in the bias values ranging within 0.02 to 0.8 mV (see Table 2). Working again on numerical simulations, we found that the bias values are highly correlated with the values of the scale factors, which makes their estimation very difficult. Our simulations showed that this correlation decreases quite quickly as the range of variation in the GraviMob orientation angles increases. More specifically, the simulations showed that the deviations of the estimated internal parameters from their reference values are less than 10 6 and the standard deviation of the residuals less than 10 6 mGal when the orientation range of GraviMob is extended to all directions. In this case, the tilt angles to which the GraviMob must be exposed during the adjustment process have to range from 180 to + 180 degrees in the three orthogonal directions of space. Another condition is that the measuring range of the accelerometers must be adjusted so that all the specific forces they experience during the experiment can be measured without entering saturation, and without adding an offset voltage. It is therefore planned to modify the measurement range and the arrangement of the accelerometers to meet these two conditions in the future version of the GraviMob system.
As seen earlier, reducing the level of residual dispersion requires multiplying the number of orientations given to the GraviMob system during the adjustment manipulation. Another fundamental point is the determination of the equations governing the variations of the internal parameters as a function of temperature. The GraviMob system is equipped with a total of seven temperature sensors. Each accelerometer has its own built-in temperature sensor and a final sensor is placed directly on the electronic acquisition board inside the glass sphere that houses it. To adjust and maintain the temperature inside the sphere during the adjustment process, it is necessary to operate with the GraviMob system installed in a temperature-controlled climate chamber. We therefore carried out a second series of laboratory measurements using the following experimental protocol:
  • installation of the GraviMob system in the climate chamber (Figure 4) in a given orientation;
  • setting a temperature set point and recording the voltages delivered by the six accelerometers measuring specific force components and internal temperatures for 5 min, after stabilisation of the internal chamber temperature;
  • repeating of step 2 after increasing the previous setpoint temperature by 0.5 °C until the entire desired temperature range has been covered;
  • repeating of the manipulation from step 1 by changing the orientation of the GraviMob until a significant number of orientations are obtained.
The first step was to calibrate the seven temperature sensors of the GraviMob system. The desired temperature range was set from 4 °C to 21 °C to include the deep water temperatures found in the Atlantic Ocean (4 °C) and the Mediterranean Sea ( 14.5 °C). At a rate of one recording every 0.5 s (2 Hz), this manipulation made it possible to obtain 600 voltages per temperature for 14 different orientations of the GraviMob system. It thus became possible to obtain 600 estimates of each of the 11 internal parameters (3 scale factors, 5 misalignment components, 3 biases) for a given temperature. The value of each parameter was then determined by averaging these 600 measurements. The standard deviation of the mean was then obtained by dividing the empirical standard deviation of the 600 measurements by the square root of 600 (≈24). The estimates of GraviMob’s internal parameters obtained by this method are presented in Figure 5, Figure 6 and Figure 7. In order to obtain mathematical functions describing the variation of these parameters with temperature, polynomial fits were performed and validated by means of chi-squared statistical tests. Our own results confirm the strong variations of the parameters with temperature, which make it essential to calibrate the system in temperature beforehand. The values of the parameters must then be set according to the mathematical model based on the working temperature of the GraviMob system measured inside the sphere.

2.3.4. Orientation of the GraviMob System Inside the Submersible

The adjustment of GraviMob’s internal parameters ends with the determination of the matrix C s o b , which fixes the orientation of the gravimob system in the submersible. The b-frame and s o -frame being orthogonal, the determination of matrix C s o b can be performed using the orthogonal Procustes method mentioned above. For this purpose, specific force measurements expressed jointly in the b-frame and s o -frame must be acquired. To this end, the adjustment manipulation consists this time in giving the GraviMob system different orientations once it is fixed in the AUV. This manipulation was carried out in the IFREMER laboratory onshore using a lifting crane that allows the AUV to be given different orientations according to the roll, pitch and yaw axes (c.f. Appendix A). By doing so, 2 N measurement pairs a M α b , a M α s o i and a M β b , a M β s o i for i = 1 , 2 , , N could be collected. During this manipulation, the INS was switched on to track the attitude angles of the AUV. We will now reason with the data for triad α , knowing that the latter is identical in every respect for triad β . The specific forces a M α s o can be obtained from the voltages V s α delivered by the accelerometers from Equation (24) as:
a M α s o = C s o s α 1 K α 1 V s α V 0 s α ,
which is possible once the scale factor matrix K α , the biais vector V 0 s α and the elements of matrix C s o s α and have been determined. The use of the INS during the adjustment process allows access to the measurements of the specific force a C b at point C, that is the AUV reference point. Using Equation (6), this specific force may be related to the gravity vector g C n at point C as:
a C b = C n b C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e + Ω ˙ i e e X P e + Ω i b b Ω i b b + Ω ˙ i b b L C b C n b g C n ,
where L C b = P C b . Similarly, the specific force a M α b at point M α is expressed as:
a M α b = C n b C e n X ¨ P e + 2 Ω i e e X ˙ P e + Ω i e e Ω i e e + Ω ˙ i e e X P e + Ω i b b Ω i b b + Ω ˙ i b b L α b C n b g M α n ,
where L α b = P M α b . Subtracting members from Equations (43) and (42) yields:
a M α b a C b = Ω i b b Ω i b b + Ω ˙ i b b L α b L C b + C n b g C n g M α n .
This equation indicates that the difference between the two specific forces depends on a term related to the lever arm between point M α and point C and the difference in the gravitational fields measured at point M α and point C respectively. The former can be bounded above using the following inequality:
Ω i b b Ω i b b + Ω ˙ i b b L α b L C b ω i b 2 + ω ˙ i b L α b L C b ,
where ω i b denotes the angular rate of the rotation of the b-frame with respect to the i-frame and ω ˙ i b its first derivative with respect to time. During static manipulations, this angular rate of rotation corresponds to that of the Earth’s rotation, and so we can pose:
ω i b = ω e = 7.292115 × 10 5 rad s 1 and ω ˙ i b = ω ˙ e 10 18 rad s 2 ,
which yields to:
Ω i b b Ω i b b + Ω ˙ i b b L α b L C b M L α b L C b with M = 0.54 μ Gal · m 1 .
The lever arm L α b L C b is known, since it was determined from a 3D digital model of the implementation of the GraviMob system and other sensors in the AUV, and can be expressed in the b-frame as:
L α b L C b = C M α b = + 2587.9 4.3 + 248.1 b ,
where the values of lever arm vector components are given in millimeters. Hence, L α b L C b = 2599.8 mm < 2.6 m . The lever arm term is therefore less than 2.6 × 0.54 = 1.4 μ Gal , and can therefore be neglected without jeopardising the targeted accuracy at 1 mGal. As for the difference in gravitational fields, its rigorous estimation presupposes knowledge of the gravity gradients in the three directions of the n-frame. If we take the value of the free air gravity gradient, i.e., 0.3 μ Gal . m 1 , as the upper limit of the amplitude of the gravity gradients, we obtain roughly
g C n g M α n < 2.6 × 0.3 = 0.8 mGal .
We therefore decided to disregard this difference even if this decision is more debatable here. If ever higher accuracy accelerometers were to be used on the GraviMob system, a mapping of the gravity gradients in the laboratory site where the adjustment manipulation is performed would become absolutely essential. By setting a M α b i a C b i for all i = 1 , 2 , , N , the orthogonal Procrustes method leads to:
C s o b = V U T ,
where U and V are the two orthogonal matrices involved in the singular value decomposition U Σ V T of the cross-correlation matrix S α defined by:
S α = i = 1 N a M α s 0 i a M α b i T .
The method was applied to data acquired for both triads in the IFREMER laboratory for 31 AUV tilts. The results given in Table 3 correspond to the Euler angles derived from the expression (A17) of matrix C s o b .

2.4. Method for Processing GraviMob-Measured Data

The gravity vector components are estimated through an Unscented Kalman Filter (UKF) [34]. This version of the Kalman filter has been preferred to the extended version (EKF) because it does not require repeating at each iteration, the determination and evaluation of the Jacobian matrix which contains the partial derivatives of the observation function in relation to the parameters sought. This filter is defined by an evolution model and an observation model. The evolution model links the parameters at a given time to the one at the previous time. By denoting X k the vector of parameters at the time k T e where T e is the sampling period, its components are given by:
X k = λ λ ˙ λ ¨ φ φ ˙ φ ¨ h h ˙ h ¨ δ δ ˙ δ ¨ χ χ ˙ χ ¨ η η ˙ η ¨ g e g ˙ e g ¨ e g n g ˙ n g ¨ n g u g ˙ u g ¨ u k T ,
where the 27 different quantities are defined as follows:
  • λ , φ , h are respectively the longitude, latitude and ellipsoidal height of point P, λ ˙ , φ ˙ , h ˙ and λ ¨ , φ ¨ , h ¨ , their first and second time derivatives;
  • δ , χ , η are respectively the yaw, pitch and roll attitude angles of the AUV, δ ˙ , χ ˙ , η ˙ and δ ¨ , χ ¨ , η ¨ , their first and second time derivatives;
  • g e , g n , g u are respectively the east, north and vertical components in the n-frame of the gravity vector g ¯ P , g ˙ e , g ˙ n , g ˙ u and g ¨ e , g ¨ n , g ¨ u , their first and second time derivatives.
The evolution model may then be expressed mathematically as:
X k + 1 = A X k + Q
where X k and X k + 1 are respectively the parameter vectors at successive times k and k + 1 , A, the evolution matrix written as a discrete Wiener stochastic process, and Q the covariance matrix that gives uncertainty to this evolution model. In Equation (49), A is a square matrix of order 27 which can be written, according to Bar-Shalom et al. [35], as:
A = E 0 9 0 9 0 9 E 0 9 0 9 0 9 E ,
where 0 9 stands for a square null matrix of order 9 and E is the square matrix of order 9 defined by:
E = 1 T e T e 2 / 2 0 0 0 0 0 0 0 1 T e 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 T e T e 2 / 2 0 0 0 0 0 0 0 1 T e 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 T e T e 2 / 2 0 0 0 0 0 0 0 1 T e 0 0 0 0 0 0 0 0 1 .
The evolution model assumes that the AUV moves at a constant speed when performing gravity surveys. Its attitude variations along the yaw, pitch and roll axes are also assumed to evolve at constant angular rotation rates. Moreover, spatial variations in gravity mean that the values of the components g e , g n , g u recorded by the GraviMob system change at each time step k due to the movement of the AUV. These components can therefore be considered as functions of time. With this in mind, the evolution model assumes that the second time derivatives g ¨ e , g ¨ n , g ¨ u are constants at all times. Given the uncertainty that affects these assumptions, the evolution model states that the following quantities:
  • λ ¨ , φ ¨ , h ¨ ,
  • δ ¨ , χ ¨ , η ¨ ,
  • g ¨ e , g ¨ n , g ¨ u ,
are all stochastic processes corresponding to white noise. It should be noted that setting all these quantities to zero at all times would result in a deterministic but unrealistic evolution model. White noise is a stochastic process with a mean of zero, and whose dispersion is fixed by its variance. Let σ λ ¨ 2 , σ φ ¨ 2 , σ h ¨ 2 be the respective variances of λ ¨ , φ ¨ , h ¨ , σ δ ¨ 2 , σ χ ¨ 2 , σ η ¨ 2 be those of δ ¨ , χ ¨ , η ¨ and σ g ¨ e 2 , σ g ¨ n 2 , σ g ¨ u 2 those of g ¨ e , g ¨ n , g ¨ u . Then, the covariance matrix Q associated with the evolution model may be expressed as [35]:
Q = F Q ( σ λ ¨ , σ φ ¨ , σ h ¨ ) 0 9 0 9 0 9 F Q ( σ δ ¨ , σ χ ¨ , σ η ¨ ) 0 9 0 9 0 9 F Q ( σ g ¨ e , σ g ¨ n , σ g ¨ u ) ,
where F Q ( σ 1 , σ 2 , σ 3 ) is the square matrix of order 9 as a function of σ 1 , σ 2 , σ 3 defined by:
F Q ( σ 1 , σ 2 , σ 3 ) = σ 1 2 T e 4 / 4 σ 1 2 T e 3 / 2 σ 1 2 T e 2 / 2 0 0 0 0 0 0 σ 1 2 T e 3 / 2 σ 1 2 T e 2 σ 1 2 T e 0 0 0 0 0 0 σ 1 2 T e 2 / 2 σ 1 2 T e σ 1 2 0 0 0 0 0 0 0 0 0 σ 2 2 T e 4 / 4 σ 2 2 T e 3 / 2 σ 2 2 T e 2 / 2 0 0 0 0 0 0 σ 2 2 T e 3 / 2 σ 2 2 T e 2 σ 2 2 T e 0 0 0 0 0 0 σ 2 2 T e 2 / 2 σ 2 2 T e σ 2 2 0 0 0 0 0 0 0 0 0 σ 3 2 T e 4 / 4 σ 3 2 T e 3 / 2 σ 3 2 T e 2 / 2 0 0 0 0 0 0 σ 3 2 T e 3 / 2 σ 3 2 T e 2 σ 3 2 T e 0 0 0 0 0 0 σ 3 2 T e 2 / 2 σ 3 2 T e σ 3 2 .
The numerical values of the standard deviations σ λ ¨ σ φ ¨ , σ h ¨ , σ δ ¨ , σ χ ¨ , σ η ¨ were estimated from a sample of position and attitude underwater data acquired during bathymetric surveys, carried out with the AUV under the same conditions as for the gravity surveys (same speed and depth). For this purpose, the position and attitude data were derived twice by a numerical method using spline functions [36]. Then the standard deviations were calculated from their statistical distributions. The estimation of the standard deviations σ g ¨ e , σ g ¨ n , σ g ¨ u was carried out using a mathematical model of the underwater gravity field, compatible with the environments traditionally explored by underwater gravimetry [37]. The values of gravity components g e , g n , g u were then calculated from the model at the positions of the same underwater data sample, and the standard deviations obtained from the same method as those derived from position and attitude data. All the values required to calculate the variances used in the evolution model are tabulated in Table 4.
The observation model relates the vector of parameters X k at the instant k to measurements from both the GraviMob system and the AUV’s position and orientation sensors. These measurements form the observation vector Z k at the instant k, which consists of the nine components defined as follows:
Z k = 1 2 a M α b + a M β b λ φ h δ χ η k T .
The choice of this observation vector implies that the position ( λ , φ , h ) and orientation ( δ , χ , η ) parameters are both observed and estimated with our data processing method. In contrast, the components of gravity are related to those of the specific forces by the non-linear relationship (20). Let us consider the function h from R 27 to R 9 defined by:
h X k = C n b C e n X ¨ P e + 2 Ω i e e X ˙ P e C n b g ¯ P n λ φ h δ χ η
The observation model can then be written as follows:
Z k = h X k + R
where R is the covariance matrix containing the uncertainties on the measurements. Assuming independent measurements, the matrix R can be expressed as:
R = σ a 2 I 3 0 0 0 0 0 0 0 σ λ 2 0 0 0 0 0 0 0 σ φ 2 0 0 0 0 0 0 0 σ h 2 0 0 0 0 0 0 0 σ δ 2 0 0 0 0 0 0 0 σ χ 2 0 0 0 0 0 0 0 σ η 2 ,
where I 3 stands for the identity matrix of order 3, σ a 2 are the variance of specific force measurements, σ λ 2 , σ φ 2 , σ h 2 are the variances of the longitude, latitude and ellipsoidal height respectively, and σ δ 2 , σ χ 2 , σ η 2 are the variances of the raw, pitch and roll angles, respectively. Given the accuracy of accelerometer measurements (1 mGal) and the uncertainties on attitude angles estimates (0.05 deg for yaw, 0.005 deg for pitch and roll; see Section 2.1), we can set the following standard deviation values:
σ a = 1 mGal ; σ δ = 0.05 deg , σ χ = σ η = 0.005 deg .
A maximum uncertainty of 2.5 m in planimetric positioning and 30 cm in altimetric positioning results in the following standard deviations in longitude, latitude and ellipsoidal height:
σ λ = 3.07 × 10 5 deg , σ φ = 2.25 × 10 5 deg , σ h = 0.30 m .
Given the matrices A, Q, R and the function h involved in the evolution model (49) and observation model (55), the application of the unscented Kalman filter becomes possible provided that we fix an initial vector X 0 and its covariance matrix. The detailed formulation of the unscented Kalman filter can be found in Simon J. Julier and Jeffrey K. Uhlmann [34]. Its application to underwater mobile gravimetry is specifically described in Roussel et al. [37].

3. Results

The performance of the GraviMob system was assessed from an underwater gravimetric survey carried out in the Mediterranean Sea in May 2016 (Figure 8b) [38]. During the survey, the AUV ran the four profiles located in the eastern and western zones shown in Figure 8a,c under different survey conditions varying depth and speed and alternating between constant depth and terrain follow-up survey modes. The surface gravimetric signal was previously determined on those profiles by the team of the French Marine Hydrographic and Oceanographic Service (SHOM). Six days of mission made it possible to complete 26 routes covering a cumulative distance of 169 km . Of those 26 routes, only 10 were found to be directly usable. The other routes require the removal of measurement periods during which the accelerometers reached saturation and/or the recording of navigation data was interrupted. In the remainder of the article, we will focus on two opposite routes taken on a profile named “910-2-S2006-076_TRANSIT” (see Figure 8b) of 9 km length. The submersible twice followed this profile, under the same navigation conditions at a constant depth of 1900 m with an average speed of 1.5 m · s 1 . From now on, the two routes will be designated by the numbers 8 and 9 respectively.
Monitoring the temperature variations inside the sphere allowed us to verify that the temperature had stabilised at the time the GraviMob measurements were taken (Figure 9). More specifically, it was observed that the internal temperature variations of the GraviMob system during the surveys on routes 8 and 9 remained below 0.04 °C and 0.004 °C respectively. With a temperature sensitivity of 15 mGal/°C, the impact on GraviMob measurements is estimated at less than 1 mGal. The hypothesis that the underwater environment acts as a thermostat keeping the temperature of the GraviMob system constant has therefore been verified.
The surface gravimetric signal provided by the SHOM consists of the free air gravity anomaly, that appears strongly smoothed compared to bathymetry (Figure 10). The weakness of the surface signal due to the high depths encountered, added to the practical low-pass filters in marine gravimetry, most probably explains the excessive smoothing of the free air anomaly towards the high frequencies. It is also worth noting the low signal variation, which does not exceed ten milligals.
The estimation of the three components of the gravity vector at the output of the UKF is represented in Figure 11. The convergence of the formal error of the UKF deduced from the covariance matrix of the estimated parameters is reached beyond 1000 computation times, which corresponds to a distance along the profile of about 750 m (Figure 12). This error represents the level of consistency of the estimated components of the gravity field given the modelling of the problem and the uncertainties in the quantities involved in the estimation process. However, it does not correspond to the accuracy of the values obtained for the gravity field components. As the formal error on the g u p component is a factor of four smaller than for the two others, the system modelling, including the evolution and observation models and their uncertainties, seems to be more appropriate for this vertical component.
Examination of the variations in gravity components (Figure 11) shows that the filtered signal still contains the contribution of non-gravitational accelerations. This finding indicates that an additional filtering operation will be necessary to attenuate these contributions. Before doing so, the mean values of components g e , g n and g e on route 8 and 9 were compared with each other. The mean values of the east and north gravity components were found to be non-identical in the two runs due to the presence of a residual bias. These components were then corrected for this bias at the output of the UKF. This correction did not change the vertical component. After correction, the mean values g e and g n of g e and g n for both routes were respectively estimated to be:
g e = 48.5 mGal g n = 38.3 mGal
The mean values g u of the vertical component g u following the two runs were respectively calculated as:
Route 8 : g u = 98 , 0160.3 mGal Route 9 : g u = 98 , 0169.3 mGal
The difference of 9 mGal observed between the two mean values can be interpreted as an estimate of the measurement precision at this stage of the processing. This value results mainly from the effect of the remaining non-gravitational accelerations on the vertical component of gravity and, to a lesser extent, from uncertainties in the determination of the internal parameters of GraviMob.
In order to compare the two signals obtained with the SHOM reference free air anomaly, a normal gravity field model at zero altitude [28] was first subtracted from the estimated gravity vector intensity g ¯ n in order to calculate a gravity anomaly δ g . A moving window filtering of length L = 3000 m allowed to smooth out the high-frequency variations affecting the two newly formed free air anomalies (Figure 13). This operation is not sufficient to make the GraviMob sea floor gravity anomalies and the SHOM reference anomalies superimposed. On average, the two were observed to be separated by about 210 mGal . Among the causes of the observed difference is the influence of the water layer above the submersible due to the fact that measurements were carried out at 1900 m depth. Correcting measurements from the gravitational effect of the water layer between underwater topographic and sea surfaces would require an accurate model of the underwater relief in order to calculate its gravimetric effect at the depth of immersion. In order to compare the sea floor and reference gravity anomalies, we opted to examine signal trends. It was therefore considered to adjust a polynomial model on each of the measured and reference gravity anomalies. A second-order polynomial model for the reference gravity anomalies was found to be appropriate. The seafloor gravity anomalies were therefore modelled by a polynomial model of the same order (Figure 14). Since the above-mentioned effects are difficult to quantify at each point of the gravity profiles, it was decided to estimate a transfer function that would best bring the trend of undersea measured gravity anomalies to the trend of reference gravity anomalies. The transfer function thus obtained makes it possible to calculate the corrections to be applied to any other gravity anomalies measured on the same profile, given similar navigation conditions. It includes both the gravity anomaly corrections needed for the upward continuation and those that compensate for the use of GraviMob’s unadapted, even misestimated internal parameters.
Considering an affine transfer function, the expression of the best transformation law that brings the gravity anomalies δ g measured during route 8 and expressed in milligals to the reference gravity may be written as:
f 8 ( δ g ) = 0.21 × δ g + 0.62
Once this function is applied on the anomalies measured from route 9, the residuals between the corrected anomalies, denoted f 8 ( δ g 9 ) , and the reference anomalies can be calculated. The histogram of these residuals (Figure 15) shows that they lie within the range between 0.5 and 3.5 mGal . The same operation can be performed using the anomalies measured during route 9 to estimate the transfer function. In that case, it may be expressed as:
f 9 ( δ g ) = 0.16 × δ g 14.8
Applied to the anomalies acquired during route 8, the residuals between the corrected anomalies, denoted f 9 ( δ g 8 ) , and the reference anomalies lie this time within the range between 0.5 and 3.5 mGal (Figure 16). These results show that, on this profile, the maximum measurement error on the estimation of gravity anomaly trends is between 3 and 4 mGal.
Equations (57) and (58) allow us to estimate the mean values of the components of the deflection from the vertical η and ξ in the east–west and north–south directions respectively by:
η = g e g and ξ = g n g ,
where g = g e 2 + g n 2 + g u 2 . Hence we obtain:
η 10 and ξ 8 .
Noting σ η , σ ξ and σ g the respective uncertainties on η , ξ and any one of the components g e or g n , then it holds that:
σ η = σ ξ σ g g ,
which in our case gives σ g g = 2 for σ g = 10 mGal . The mean values of the deflection of vertical were compared to those calculated from recent geopotential models (cf. Table 5) on the studied profile. For this purpose, the mean values of east–west and north–south deflection components at 1000 points evenly distributed along the studied profile were determined by different models (cf. Table 6). The geopotential model values were calculated from the International Centre for Global Earth Models (ICGEM) [40].
The results obtained (Table 6) show the good adequacy of the value of the north–south component ξ (-8 arcsec) with that given by the geopotential models. However, the east–west component η was clearly overestimated (10 arcsec compared to 3–4 arcsec), which would correspond to an error on the vertical gravity component of 30 mGal to 35 mGal (cf. Equation (63)). This error is incompatible with the results which indicate agreement between the vertical gravity component values to within 9 mGal (cf. Equation (58)). As mentioned before, the gravity components g e and g n were corrected for a residual bias. This bias can be represented by a disruptive horizontal acceleration whose direction is close to that of the AUV trajectory along the studied profile. Given its direction, the profile studied is only inclined 36° to the east–west direction. Therefore, the bias will affect the estimation of the gravity component g e more strongly. More precisely, a perturbing horizontal acceleration of 40 mGal along the profile could explain a deviation of 32 mGal, that 40 / cos 36 , on the east gravity component. This result shows that the calculation of the bias correction from only the measurements obtained on a direct and reverse run along the studied profile is not yet sufficiently accurate, and must therefore be refined by increasing the number of runs and profiles. The bias is likely to be due to an erroneous estimate of the orientation of the sensitive axes of the accelerometers relative to those of the vehicle frame. It can be reduced or even eliminated by refining the estimation of GraviMob’s internal parameters.

4. Discussion

The comparison method used here to quantify the accuracy of the GraviMob system appears to be the most rigorous given the reference data available to us. The SHOM surface reference gravity data is certainly not ideal but it is the only existing quality gravity data that can be used in this case. A bathymetric model of the metric resolution, not available in this experiment, would require significant work of bathymetric survey but would nevertheless facilitate the upward continuation of observed data. That said, the repeatability of the gravimetric signal measurement gives encouraging results, which are likely to improve significantly with changes to the instrument design and data processing method. We therefore propose to review the changes that we consider to be the most important for improving the performance of the GraviMob system.
As seen in Section 2.3 and Section 3, the reliable estimation of the internal parameters of the GraviMob system and their temperature dependence is a crucial point for the processing of the data acquired by the system. This determination requires a sufficiently wide range of permitted measurements so that the GraviMob system can be oriented in many different directions during the parameter adjustment phases without the accelerometers becoming saturated. This can be achieved by using pyramid-shaped triads (Figure 17), such that the inclinations of the axes of the three accelerometers with respect to the vertical direction are identical when the instrument is at rest. By doing so, a measurement range of ± g with a sensitivity of 0.2 mGal can be achieved without having to add an offset voltage to any accelerometer. A sensitivity of 0.2 mGal instead of 0.02 mGal is sufficient due to the uncertainty in the accelerometer measurements. This modification is essential to reduce, if not cancel, the residual bias that has strongly affected the estimation of the horizontal gravity components by our calculation method.
As seen in Section 2.2, removing the lever arm that separates the gravity measurement point from the position measurement point would greatly simplify and make the processing method more reliable. This would require an assembly that would allow an INS and the GraviMob system to be arranged concentrically. This solution is possible through a radical change to the current design of the instrument, which will necessarily involve the manufacturers of inertial navigation systems.
For the data processing method based on the unscented Kalman filter, a first way to reduce the error bar on the gravity components is to use more accurate position and orientation measurements. A second way to reduce the uncertainty of the gravity components is to modify the evolution model used in the Kalman filter to make it more physically realistic. As seen in Section 2.4, the evolution model used in the Kalman filter assumes that the second derivatives of the gravity components are constant within one white noise process. This assumption amounts to considering that the covariance function of gravity would have an unbounded support. However, it is known that the covariance function of real gravity is a decreasing function whose rate of decrease depends on a finite correlation length. One solution to this problem is to adjust the parameters of the stochastic Markov model that describes the gravity components to ensure that the corresponding covariance function coincides with an empirical covariance function determined for the study area. The method for performing such an adjustment has been described theoretically and simulated numerically in Verdun et al. [45], based on an idea presented in the book by Jekeli [24]. The use of this new evolution model should lead to a substantial improvement in the estimation of the gravity components from the same position and orientation measurements as in our calculation.

5. Conclusions

In conclusion, we have presented the GraviMob system, an innovative and lightweight inertial gravity sensor that can be easily installed on an AUV to perform undersea gravity surveys. GraviMob is a vector gravimeter that allows the determination of the three components, east, north and vertical of gravity. The sensor is formed by two parallelepiped triads each containing three electrostatic accelerometers whose geometric axes are mutually orthogonal. The sensor is enclosed in a sealed spherical enclosure which is then installed in the AUV head. A data acquisition system associated with a digitiser is used to record the measurements of the six accelerometers and the internal temperature measurements.
The measurement principle detailed in the article is based on the combined processing of the GraviMob system measurements with the position and orientation measurements provided by the AUV’s own navigation system. The measurement processing involves the prior determination of internal parameters of the GraviMob system, including the misalignments of the sensitive axes of the accelerometers, the accelerometer biases and scale factors, and the orientation of the accelerometer axes with respect to those of the AUV’s inertial navigation system. A comprehensive experimental method for estimating the values of the internal parameters with high accuracy, described in detail in the article, has been developed. A calibration of the internal parameters according to temperature was also carried out. Further processing of the measurements is based on the use of an unscented Kalman filter which does not require linearisation of the GraviMob’s observation equation.
A trial cruise in the Mediterranean sea tested the performance of the GraviMob system on board an AUV operated by IFREMER immersed at a depth of 1900 m. The validation of GraviMob’s measurements consisted in comparing them with reference measurements, namely, surface gravity measurements carried out by the SHOM prolongated downwards and gravity values calculated from recent global gravity geopotential field models. The results of an analysis carried out on a profile run in two opposite directions indicate that the deviations of the vertical component of gravity from the reference remain below 4 mGal. The east and north components of the vertical deflection deduced from the GraviMob system measurements were compared with those computed from recent global gravity field models at the AUV measurement depth. Although the north component of deflection agrees remarkably well with that given by the models ( 8 ± 2 arcsec), there is a significant discrepancy on the east component ( 10 ± 2 arcsec compared to 3 to 3.5 arcsec on average). We believe that this discrepancy can be compensated for by improving the estimation of the orientation angles of the GraviMob’s sensitive axes with respect to those of the AUV’s inertial navigation system. The replacement of parallelepiped triads with pyramid triads in the GraviMob system is the key to achieving this goal.
In addition to changing the shape of the triads, further development of the GraviMob system will involve improving the evolution model used in the Kalman filter to take into account the natural covariances of the gravity field. The performance of the GraviMob prototype is undoubtedly very promising and could be further improved by using accelerometers with better resolution or by enhancing the AUV pose estimation. With the planned improvements, we hope to obtain an estimate of the three components of gravity with an uncertainty of less than 2 mGal.

Author Contributions

Conceptualisation of the project: M.M., J.V. and J.C.; funding acquisition: M.M. and J.V; project administration: M.M.; head of sea trial missions: M.M. and C.P.; participants in sea trials: M.M., C.P., C.R., J.V., O.K. and J.C.; supervision of C.R.: J.V., J.C. and M.M.; supervision of O.K.: J.V., J.C., C.P. and J.A.; technical guidance to software implementation and computer management: F.D.; algorithm design and software implementation: C.R., O.K. and J.V.; gravity data processing: C.R., J.V. and J.C.; navigation data processing: M.-É.B.; design of the sensor: J.C., J.V. and J.-F.D.; assembly and integration of the sensor: J.-F.D.; design and implementation of the data acquisiton software: J.-F.D.; design and implementation of experimental protocols for the estimation of the internal parameters of the system: C.R., O.K., J.C., J.V., J.A. and C.P.; GraviMob’s internal parameter computation: C.R., O.K., J.V. and J.C.; result analysis and discussion: J.V., J.C., C.R., O.K., C.P., J.A. and M.M.; writing and editing: J.V., C.R., M.M. and J.C. All authors have read and agreed to the published version of the manuscript.

Funding

The different stages of the GraviMob project were successively financed by CNRS INSU-CESSUR, CPER Région Bretagne and Labex MER. The sea trials were financed by the French Oceanographic Fleet (FOF). Fundings were also obtained from the laboratories involved in the project (Ocean Geosciences Lab from CNRS, UBO, UBS and GeF Lab from the Cnam). Clément Roussel’s thesis was financed by the French Ministry of Defence (DGA) and the Pays de la Loire Region. Cnam funded the internship of Ossama Kharbou.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Maps were obtained using GMT software (https://www.generic-mapping-tools.org/ accessed on 12 March 2022). The geopotential models used in the paper were obtained from the ICGEM (http://icgem.gfz-potsdam.de/home accessed on 12 March 2022).

Acknowledgments

We would like to thank the engineers and technical staff of IFREMER Toulon systèmes sous-maris for the professionalism and thoroughness with which they conducted the GraviMob marine cruise. We are also very grateful to Captain F. Lofficial and his crew for welcoming us so warmly on board the vessel Europe. We are grateful to anonymous reviewers who helped improve the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
AUVAutonomous Underwater Vehicle
EKFExtended Kalman Filter
FOFFrench Oceanographic Fleet
GRACEGravity Recovery and Climate Experiment
GraviMobMobile Gravimetry System
GNSSGlobal Navigation Satellite Systems
GOCEGravity field and steady-state Ocean Circulation Explorer
ICGEMInternational Centre for Global Earth Models
IFREMERFrench National Institute for Ocean Science
IGNFrench National Institute of Geographic and Forest Information
INSInertial Navigation System
JAMSTECJapan Agency for Marine-Earth Science and Technology
LAGEOSLAser GEOdynamics Satellite
LIMO-gLIght MObile gravimetry system
LODLength Of Day
ROVRemotely Operated Vehicle
UKFUnscented Kalman Filter
USBLUltra-Short Baseline Position System
UUVUnmanned Underwater Vehicle
SHOMFrench Marine Hydrographic and Oceanographic Service

Appendix A. Definition of Useful Coordinate Reference Frames

In the paper, use is made of six coordinate reference systems designated respectively by i-frame, e-frame, n-frame, b-frame, s-frame and s o -frame. These are defined as follows.
  • The i-frame (inertial frame) is a non-rotating, inertial reference frame which has its origin at the centre O of the Earth and axes aligned with the directions of fixed stars (Figure A1). This frame consists of an orthogonal, right-handed axis set defined by the axes O x i , O y i and O z i with O z i coincident with the Earth’s polar axis assumed to be invariant in direction. The i-frame is the only reference frame where Newton’s Second Law of motion is valid.
  • The e-frame (Earth frame) is the Earth fixed, reference frame used for location definition. Its origin is at the centre O of the Earth and axes are fixed with respect to the Earth (Figure A1). This frame consists of an orthogonal, right-handed axis set O x e , O y e , O z e typically defined with O z e parallel to the Earth’s polar axis and O x e lying along the intersection of the Greenwich meridian with the Earth’s equatorial plane. The e-frame rotates with respect to the i-frame at a angular rate ω e about the axis O z e .
  • The n-frame (navigation frame) is a local geographic frame which has its origin C at the point whose position is effectively measured by the AUV navigation system (Figure A1). The frame consists of an orthogonal, right-handed axis set whose axes C x n , C y n and C z n are aligned respectively with the directions of east, north and the upward normal to a reference ellipsoid passing through the point C. Thus, the n-frame moves with the submersible and the movement of its origin C can be determined with respect to the e-frame using the AUV navigation system. The components of the gravity vector g ¯ n are usually expressed in the n-frame, thus giving respectively its east g e , north g n and vertical g u components.
Figure A1. View of the inertial reference frame O x i y i z i (i-frame), the Earth fixed, reference frame ( O x e y e z e ) (e-frame) and the navigation reference frame ( C x e y e z e ) (n-frame). The point O corresponds to the centre of the Earth represented by an ellipsoid of revolution whose axis of revolution is the pole axis ( O z i ) . The number triplet ( λ , φ , h ) denote respectively the longitude, the geographic latitude and the ellipsoidal height of the point C.
Figure A1. View of the inertial reference frame O x i y i z i (i-frame), the Earth fixed, reference frame ( O x e y e z e ) (e-frame) and the navigation reference frame ( C x e y e z e ) (n-frame). The point O corresponds to the centre of the Earth represented by an ellipsoid of revolution whose axis of revolution is the pole axis ( O z i ) . The number triplet ( λ , φ , h ) denote respectively the longitude, the geographic latitude and the ellipsoidal height of the point C.
Remotesensing 14 02513 g0a1
4.
The b-frame (body frame) consists of a orthogonal, right-handed axis set which axes C x b , C y b and C z b are aligned respectively with the pitch, roll and yaw axes of the submersible. Its orientation with respect to the n-frame is used for defining the attitude of the submersible (Figure A2). In the paper, we assume that the b-frame and the n-frame have both the same origin C. The point C will henceforth be called the AUV reference point.
5.
The s-frame (sensor frame) is an acceleration sensor coordinate frame with axes parallel to the sensor input axes of one given accelerometer (Figure A2). These axes are non-coplanar, but not necessarily orthogonal depending on both the design of the accelerometer supporting triad and the misalignments which affect unavoidably the sensor input axis directions. Each triad of accelerometers defines its own s-frame the origin of which is at the intersection of its three input axes. In the paper, we denote by M α (resp. M β ) the origin of the s-frame s α (resp. s β ) defined by the triad of accelerometers labelled by α (resp. β ).
6.
The s o -frame (orthogonal sensor frame) is also an acceleration sensor coordinate frame whose origin is at M α (resp. M β ) for the triad α (resp. β ). Unlike the s-frame, the s o -frame is an orthogonal, right-handed reference frame whose x-axis is parallel to the x-axis of the s-frame (Figure A2). To be more specific, let e x s α , e y s α , e z s α be a set of three unit vectors of the s α -frame axes x s α , y s α and z s α respectively. Let e x s o , e y s o , e z s o be a set of three orthogonal unit vectors of the s o -frame axes x s o , y s o and z s o respectively. The vectors e x s α o , e y s α o , e z s α o can be defined as follows:
  • e x s α o = e x s α ;
  • e y s α o belongs to the plane defined by the two vectors e x s α and e y s α whilst being orthogonal to the vector e x s α o ;
  • e z s α o = e x s α o × e y s α o .
The same reasoning can be used to define the frame s β for the second triad. The s o -frame is essentially used within the processing to determine GraviMob’s internal parameters (cf. Section 2.3).
Figure A2. (a) View of triplets of unit vectors (in red) which defined respectively the sensor reference frames for the triads α (top) e x s α , e y s α , e z s α and β (bottom) e x s β , e y s β , e z s β . The shaded areas indicate the arrangement of the accelerometers in the two parallelepiped triads. The misalignment angles have been deliberately exaggerated to make the figure easier to read. The corresponding triplet of unit, mutually orthogonal vectors (in black) which defined the orthogonal sensor frames for the triads α (top) and β (bottom) are also represented. (b) View of relative positions of points C (AUV reference point), P (GraviMob’s sensitive point), M α and M β (centres of α and β triads respectively) in the AUV. The locations of s α , s α o , s β , s β o s-frame and b-frame reference frames are also given. ( C P ) is the revolution axis of the AUV (roll axis). The angles η , χ , δ are repectively the roll, pitch and yaw angles. It should be noted that the vectors C P , P M α and P M β can be considered as lever arms. (c) Location of the b-frame ( C x b y b z b ) in the AUV. The revolution axis ( C P ) can also be seen.
Figure A2. (a) View of triplets of unit vectors (in red) which defined respectively the sensor reference frames for the triads α (top) e x s α , e y s α , e z s α and β (bottom) e x s β , e y s β , e z s β . The shaded areas indicate the arrangement of the accelerometers in the two parallelepiped triads. The misalignment angles have been deliberately exaggerated to make the figure easier to read. The corresponding triplet of unit, mutually orthogonal vectors (in black) which defined the orthogonal sensor frames for the triads α (top) and β (bottom) are also represented. (b) View of relative positions of points C (AUV reference point), P (GraviMob’s sensitive point), M α and M β (centres of α and β triads respectively) in the AUV. The locations of s α , s α o , s β , s β o s-frame and b-frame reference frames are also given. ( C P ) is the revolution axis of the AUV (roll axis). The angles η , χ , δ are repectively the roll, pitch and yaw angles. It should be noted that the vectors C P , P M α and P M β can be considered as lever arms. (c) Location of the b-frame ( C x b y b z b ) in the AUV. The revolution axis ( C P ) can also be seen.
Remotesensing 14 02513 g0a2

Appendix B. Coordinate and Matrix Transformations

The five coordinate frames i-frame, e-frame, n-frame, b-frame and s o -frame are assumed to be defined by three consecutively numbered (or lettered) vectors which are mutually perpendicular to one another in the right-hand sense. In this case, the transformation matrix from a 2 coordinate frame to a 1 coordinate frame is defined as a square 3 × 3 matrix the columns of which are an orthogonal set of unit vectors, each equal to a unit vector along axis of coordinate frame a 2 as projected onto the axes of coordinate frame a 1 ; thus,
C a 2 a 1 = u 1 a 2 a 1 u 2 a 2 a 1 u 3 a 2 a 1 ,
where u i a 2 a 1 is the unit vector along a 2 frame axis i projected on coordinate frame a 1 axes. From this definition, it results that the element in row i, column j of C a 2 a 1 equals the cosine of the angle between frame a 1 axis i and frame a 2 axis j, and that the transpose of C a 2 a 1 equals its inverse, that is
X a 1 = C a 2 a 1 X a 2 X a 2 = C a 1 a 2 X a 1 ,
where
C a 1 a 2 = ( C a 2 a 1 ) 1 = ( C a 2 a 1 ) T ,
or, equivalently
C a 1 a 2 C a 2 a 1 = C a 2 a 1 C a 1 a 2 = I 3 ,
where I 3 is the identity matrix of order 3. Matrix C a 2 a 1 may be called more specifically the direction cosine matrix. If one of the coordinate frames is not orthogonal such as s-frame, the relation (A2) is still valid but the corresponding set of vectors is no longer orthogonal. Thus, the inverse of matrix C a 2 a 1 has to be explicitly calculated.
Let us now consider three coordinate frames a 1 , a 2 , a 3 . In that case the transformation matrix C a 1 a 3 from a 1 to a 3 may be expressed as a function of the transformation matrices from a 1 to a 2 and a 2 to a 3 respectively as:
C a 1 a 3 = C a 2 a 3 C a 1 a 2 .
The previous equation can be generalised for an arbitrary number p of reference frames a 1 , a 2 , a 3 , ..., a p 1 , a p by multiplying the successive transformation matrices as follows:
C a 1 a p = C a p 1 a p C a p 2 a p 1 C a 2 a 3 C a 1 a 2 .
The time rate change of a direction cosine matrix is important to later mathematical developments. Let C a 2 a 1 be a direction cosine matrix, and ω a 1 a 2 the angular velocity of the a 2 -frame with respect to the a 1 -frame. Then, the time rate of change C ˙ a 2 a 1 of C a 2 a 1 is given by
C ˙ a 2 a 1 = ω a 1 a 2 × C a 2 a 1 = Ω a 1 a 2 a 1 C a 2 a 1 ,
where Ω a 1 a 2 a 1 is the skew symmetric matrix associated with the angular velocity ω a 1 a 2 coordinatised in the a 1 -frame.
The skew symmetric matrix coordinatised in the a 1 -frame can be related to the skew symmetric matrix coordinatised in the a 2 -frame by using
Ω a 1 a 2 a 1 = C a 2 a 1 Ω a 1 a 2 a 2 C a 1 a 2 .
Thus, an equivalent form of Equation (A7) may be expressed as
C ˙ a 2 a 1 = C a 2 a 1 Ω a 1 a 2 a 2 .
It should be noticed that Equation (A8) is still valid for any matrices Γ a 1 and Γ a 2 expressed respectively in the a 1 -frame and a 2 -frame thus giving:
Γ a 1 = C a 2 a 1 Γ a 2 C a 1 a 2 .
The first and second-order time derivative of equation X a 1 = C a 2 a 1 X a 2 can be calculated using Equation (A9); it thus follows respectively:
X ˙ a 1 = C a 2 a 1 X ˙ a 2 + C a 2 a 1 Ω a 1 a 2 a 2 X a 2 ;
X ¨ a 1 = C a 2 a 1 X ¨ a 2 + 2 C a 2 a 1 Ω a 1 a 2 a 2 X ˙ a 2 + C a 2 a 1 Ω a 1 a 2 a 2 Ω a 1 a 2 a 2 + Ω ˙ a 1 a 2 a 2 X a 2 .

Appendix C. Explicit Forms of Transformation and Skew Symmetric Matrices

Appendix C.1. Earth Frame to Navigation Frame

By denoting λ the longitude, and φ the latitude of the AUV reference point, the transformation matrix C e n from the e-frame to the n-frame may be expressed as:
C e n = sin λ cos λ 0 cos λ sin φ sin λ sin φ cos φ cos λ cos φ sin λ cos φ sin φ .

Appendix C.2. Body Frame to Navigation Frame

The attitude of the AUV corresponds to the relative orientation of the b-frame with respect to the n-frame. This orientation is defined by the three angles of yaw ( δ ), pitch ( χ ) and roll ( η ). These angles are such that the b-frame can be deduced from the n-frame by the following successive transformations:
  • π 2 δ angle, clockwise rotation around the inital z axis;
  • χ angle, counter-clockwise rotation around the intermediate y axis;
  • η angle, counter-clockwise rotation around the final x axis.
When these three angles are zero, the axes of the b-frame are directed north, west and vertical respectively. The transformation matrix C b n from the b-frame to the n-frame is given by:
C b n = sin δ cos χ sin δ sin χ sin η cos δ cos η sin δ sin χ cos η + cos δ sin η cos δ cos χ cos δ sin χ sin η + sin δ cos η cos δ sin χ cos η sin δ sin η sin χ cos χ sin η cos χ cos η .

Appendix C.3. Orthogonal Sensor Frame to Sensor Frame

The transformation matrix C s o s α (resp. C s o s β ) from the s o -frame to the s α -frame (resp. s β -frame) can be expressed according to the definition of the s o -frame given in Appendix A. Columns of the matrix C s o s α (resp. C s o s β ) are each formed by the components of the unitary vectors along s o -frame axes expressed relative to the s α -frame (resp. s β -frame). We therefore have:
C s o s α = 1 τ x y τ x z 0 τ y y τ y z 0 0 τ z z .
The misalignment of the axes of s-frame being small, the components of the column vectors are almost those of orthogonal unit vectors. We therefore have at first order:
τ y y 1 , τ z z 1 , τ x y 0 , τ x z 0 , τ y z 0 .
The inverse C s α s o of transformation matrix C s o s α is given by Equation (37).

Appendix C.4. Orthogonal Sensor Frame to Body Frame

The transformation matrix C s o b from the s o -frame to the b-frame is essentially an orthogonal matrix that defines the relative orientation of two orthogonal right-handed reference frames. Noting θ x , θ y , θ z , the three rotation angles, also known as Euler angles, such that the s o -frame can be deduced from the b-frame by the following successive transformations:
  • θ z angle, counter-clockwise rotation around the inital z-axis;
  • θ y angle, counter-clockwise rotation around the intermediate y-axis;
  • θ x angle, counter-clockwise rotation around the final x-axis.
We then have
C s o b = cos θ z cos θ y sin θ z cos θ x + cos θ z sin θ y sin θ x sin θ z sin θ x + cos θ z sin θ y cos θ x sin θ z cos θ y cos θ z cos θ x + sin θ z sin θ y sin θ x cos θ z sin θ x + sin θ z sin θ y cos θ x sin θ y cos θ y sin θ x cos θ y cos θ x .

Appendix C.5. Earth Frame with Respect to Inertial Frame

The skew symmetric matrix Ω i e e and its first time derivative Ω ˙ i e e are respectively given by:
Ω i e e = 0 ω e 0 ω e 0 0 0 0 0 , Ω ˙ i e e = 0 ω ˙ e 0 ω ˙ e 0 0 0 0 0 ,
where ω e is the Earth angular velocity and ω ˙ e its first time derivative.

Appendix C.6. Body Frame with Respect to Inertial Frame

The skew symmetric matrix Ω i b b is associated with the angular velocity vector ω i b b . The components of this vector can be directly provided by the measures of three gyros the sensitive axes of which are mutually orthogonal. The same can also be obtained from position and orientation measurements such as those made on GraviMob from the following equation:
ω i b b = C n b C e n ω i e e + ω e n n + ω n b b ,
where, denoting by λ the longitude, φ , the latitude, δ , χ , η , the yaw, pitch and roll angles,
ω i e e = 0 0 ω e , ω e n n = φ ˙ λ ˙ cos φ λ ˙ sin φ , ω n b b = δ ˙ sin χ + η ˙ δ ˙ cos χ sin η + χ ˙ cos η δ ˙ cos χ cos η χ ˙ sin η .
Once the components ω i b , x b , ω i b , y b , ω i b , z b of ω i b b have been determined, the skew symmetric matrix Ω i b b may be expressed as:
Ω i b b = 0 ω i b , z b ω i b , y b ω i b , z b 0 ω i b , x b ω i b , y b ω i b , x b 0 .
The first time derivative Ω ˙ i b b of the matrix Ω i b b can either be calculated analytically or numerically. It should be noted that the use of two triads in the GraviMob system avoids the need to calculate the matrices Ω i b b and Ω ˙ i b b in order to process the data.

Appendix D. Constrained Optimisation

As seen in Section 2.3.2, estimating the direction cosines c x , i , c y , i and c z , i for each inclination consists in minimizing N independent cost functions f i as (34) under the constraints stating that c x , i 2 + c y , i 2 + c z , i 2 = 1 for all i varying from 1 to N. If we omit for simplicity index i, then we can formulate the minimisation problem as the determination of:
min c f ( c ) with h ( c ) = 0 ,
where c = c x c y c z and h ( c ) = c x 2 + c y 2 + c z 2 1 .
The solution proposed by Yang et al. [32] is based on Karush–Kuhn–Tucker condition [46] defined by the following equation:
L ( c , u ) = f ( c ) + u h ( c ) = 0 ,
where L is the lagrangian function of the problem and u the Lagrange multiplier associated with the constraint h ( c ) = 0 .
Equation (A21) can be iteratively solved by using the Newton method. Given an initial solution ( c 0 , u 0 ) , then we have at iteration k:
c k + 1 = c k + δ c and u k + 1 = u k + δ u ,
with ( δ c , δ u ) that satisfy the following linear system:
2 L ( c , u ) h ( c ) h ( c ) T 0 δ c δ u = L ( c , u ) h ( c ) T .
The explicit content of the linear system (A23) is obtained from the equations given below. Denoting
μ = μ x μ y μ z = g M n k x x c x + g M n k x y c y + g M n k x z c z + v x 0 g M n k y y c y + g M n k y z c z + v y 0 g M n k z z c z + v z 0 ,
the different elements of the matrices involved in Equation (A23) are expressed as follows:
h ( c ) = c x 2 + c y 2 + c z 2 1 ; h ( c ) = 2 c x c y c z ; L ( c , u ) = L / c x L / c y L / c z ; = 2 g M n k x x ( μ x v x ) + 2 u c x 2 g M n k x y ( μ x v x ) + 2 g M n k y y ( μ y y y ) + 2 u c y 2 g M n k x z ( μ x v x ) + 2 g M n k y z ( μ y v y ) + 2 g k z z ( μ z v z ) + 2 u c z ; 2 L ( c , u ) = l x x l x y l x z l y x l y y l y z l z x l z y l z z ,
with:
l x x = 2 L c x 2 = 2 g M n 2 k x x 2 + 2 u , l x y = c y L c x = 2 g M n 2 k x x k x y , l x z = c z L c x = 2 g M n 2 k x x k x z , l y x = c x L c y = 2 g M n 2 k x y k x x , l y y = 2 L c y 2 = 2 g M n 2 k x y 2 + 2 g M n 2 k y y 2 + 2 u , l y z = c z L c y = 2 g M n 2 k x y k x z + 2 g M n 2 k y y k y z , l z x = c x L c z = 2 g M n 2 k x z k x x , l z y = c y L c z = 2 g M n 2 k x z k x y + 2 g M n 2 k y z k y y , l z z = 2 L c z 2 = 2 g M n 2 k x y 2 + 2 g M n 2 k y z 2 + 2 g M n 2 k z z 2 + 2 u .

References

  1. Herzig, P.M.; Hannington, M.D. Polymetallic Massive Sulfides and Gold Mineralization at Mid-Ocean Ridges and in Subduction-Related Environments; CRC Press: Boca Raton, FL, USA, 2000; pp. 347–368. [Google Scholar]
  2. Rona, P.A. Resources of the Sea Floor. Science 2003, 299, 673–674. [Google Scholar] [CrossRef] [PubMed]
  3. Hannington, M.; Jamieson, J.; Monecke, T.; Petersen, S.; Beaulieu, S. The abundance of seafloor massive sulfide deposits. Geology 2011, 39, 1155–1158. [Google Scholar] [CrossRef]
  4. Escartín, J.; Lin, J. Ridge offsets, normal faulting, and gravity anomalies of slow spreading ridges. J. Geophys. Res. Solid Earth 1995, 100, 6163–6177. [Google Scholar] [CrossRef] [Green Version]
  5. Maia, M.; Goslin, J.; Gente, P. Evolution of the accretion processes along the Mid-Atlantic Ridge north of the Azores since 5.5 Ma: An insight into the interactions between the ridge and the plume. Geochem. Geophys. Geosyst. 2007, 8. [Google Scholar] [CrossRef] [Green Version]
  6. Maia, M.; Sichel, S.; Briais, A.; Brunelli, D.; Ligi, M.; Ferreira, N.; Campos, T.; Mougel, B.; Brehme, I.; Hémond, C.; et al. Extreme mantle uplift and exhumation along a transpressive transform fault. Nat. Geosci. 2016, 9, 619–623. [Google Scholar] [CrossRef] [Green Version]
  7. Ballu, V.; Dubois, J.; Deplus, C.; Diament, M.; Bonvalot, S. Crustal structure of the Mid-Atlantic Ridge south of the Kane Fracture Zone from seafloor and sea surface gravity data. J. Geophys. Res. Solid Earth 1998, 103, 2615–2631. [Google Scholar] [CrossRef] [Green Version]
  8. Evans, R.L. A seafloor gravity profile across the TAG Hydrothermal Mound. Geophys. Res. Lett. 1996, 23, 3447–3450. [Google Scholar] [CrossRef]
  9. Araya, A.; Kanazawa, T.; Shinohara, M.; Yamada, T.; Fujimoto, H.; Iizasa, K.; Ishihara, T. A gravity gradiometer to search for submarine ore deposits. In Proceedings of the 2011 IEEE Symposium on Underwater Technology and Workshop on Scientific Use of Submarine Cables and Related Technologies, Tokyo, Japan, 5–8 April 2011; pp. 1–3. [Google Scholar] [CrossRef]
  10. Fujimoto, H.; Kanazawa, T.; Shinohara, M.; Araya, A.; Yamada, T.; Mochizuki, K.; Ishihara, T.; Iizasa, K. Development of a hybrid gravimeter system onboard an underwater vehicle. In Proceedings of the 2011 IEEE Symposium on Underwater Technology and Workshop on Scientific Use of Submarine Cables and Related Technologies, Tokyo, Japan, 5–8 April 2011; pp. 1–3. [Google Scholar] [CrossRef]
  11. Araya, A.; Kanazawa, T.; Shinohara, M.; Yamada, T.; Fujimoto, H.; Iizasa, K.; Ishihara, T. Gravity gradiometer implemented in AUV for detection of seafloor massive sulfides. In Proceedings of the 2012 Oceans, Hampton Roads, VA, USA, 14–19 October 2012; pp. 1–4. [Google Scholar] [CrossRef]
  12. Shinohara, M.; Yamada, T.; Kanazawa, T.; Uehira, K.; Fujimoto, H.; Ishihara, T.; Araya, A.; Iizasa, K.; Tsukioka, S. Development of an underwater gravimeter and the first observation by using autonomous underwater vehicle. In Proceedings of the 2013 IEEE International Underwater Technology Symposium (UT), Tokyo, Japan, 5–8 March 2013; pp. 1–6. [Google Scholar] [CrossRef]
  13. Araya, A.; Shinohara, M.; Kanazawa, T.; Fujimoto, H.; Yamada, T.; Ishihara, T.; Iizasa, K.; Tsukioka, S. Development and demonstration of a gravity gradiometer onboard an autonomous underwater vehicle for detecting massive subseafloor deposits. Ocean Eng. 2015, 105, 64–71. [Google Scholar] [CrossRef]
  14. Shinohara, M.; Yamada, T.; Ishihara, T.; Araya, A.; Kanazawa, T.; Fujimoto, H.; Uehira, K.; Tsukioka, S.; Omika, S.; Iizasa, K. Development of an underwater gravity measurement system using autonomous underwater vehicle for exploration of seafloor deposits. In Proceedings of the OCEANS 2015, Genova, Italy, 18–21 May 2015; pp. 1–7. [Google Scholar] [CrossRef]
  15. Shinohara, M.; Kanazawa, T.; Fujimoto, H.; Ishihara, T.; Yamada, T.; Araya, A.; Tsukioka, S.; Omika, S.; Yoshiume, T.; Mochizuki, M.; et al. Development of a High-Resolution Underwater Gravity Measurement System Installed on an Autonomous Underwater Vehicle. IEEE Geosci. Remote Sens. Lett. 2018, 15, 1937–1941. [Google Scholar] [CrossRef]
  16. Ishihara, T.; Shinohara, M.; Fujimoto, H.; Kanazawa, T.; Araya, A.; Yamada, T.; Iizasa, K.; Tsukioka, S.; Omika, S.; Yoshiume, T.; et al. High-resolution gravity measurement aboard an autonomous underwater vehicle. Geophysics 2018, 83, G119–G135. [Google Scholar] [CrossRef]
  17. Kinsey, J.C.; Tivey, M.A.; Yoerger, D.R. Toward high-spatial resolution gravity surveying of the mid-ocean ridges with autonomous underwater vehicles. In Proceedings of the OCEANS 2008, Quebec City, QC, Canada, 15–18 September 2008; pp. 1–10. [Google Scholar] [CrossRef]
  18. Izraelevitz, J. Optimal trajectory generation for draped AUV gravity surveys. In Proceedings of the OCEANS 2011, Santander, Spain, 6–9 June 2011; pp. 1–8. [Google Scholar] [CrossRef]
  19. Kinsey, J.C.; Tivey, M.A.; Yoerger, D.R. Dynamics and navigation of autonomous underwater vehicles for submarine gravity surveying. Geophysics 2013, 78, G55–G68. [Google Scholar] [CrossRef]
  20. Zhang, Z.; Li, J.; Zhang, K.; Yu, R. Experimental Study on Underwater Moving Gravity Measurement by Using Strapdown Gravimeter Based on AUV Platform. Mar. Geod. 2021, 44, 108–135. [Google Scholar] [CrossRef]
  21. Saint-Jean, B.d. Étude et Développement d’un Système de Gravimétrie Mobile. Ph.D. Thesis, Observatoire de Paris, Paris, France, 2008. [Google Scholar]
  22. Roussel, C. Expérimentation d’un Gravimètre Mobile Léger et Novateur pour la Mesure du Champ de Gravité en Fond de Mer. Ph.D. Thesis, Conservatoire National des Arts et Métiers (Cnam), Paris, France, 2017. [Google Scholar]
  23. El-Sheimy, N.; Hou, H.; Niu, X. Analysis and Modeling of Inertial Sensors Using Allan Variance. IEEE Trans. Instrum. Meas. 2008, 57, 140–149. [Google Scholar] [CrossRef]
  24. Jekeli, C. Inertial Navigation Systems with Geodetic Applications; Walter de Gruyter: Berlin, Germany, 2012. [Google Scholar] [CrossRef]
  25. Rogers, R.M. Applied Mathematics in Integrated Navigation Systems, 3rd ed.; American Institute of Aeronautics & Astronautics: Reston, VA, USA, 2007. [Google Scholar]
  26. Titterton, D.; Weston, J.; Weston, J. Strapdown Inertial Navigation Technology; IEE Radar, Sonar, Navigation and Avionics Series; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2004. [Google Scholar]
  27. Cai, S.; Zhang, K.; Wu, M.; Huang, Y. Long-Term Stability of the SGA-WZ Strapdown Airborne Gravimeter. Sensors 2012, 12, 11091–11099. [Google Scholar] [CrossRef]
  28. Moritz, H. Geodetic Reference System 1980. J. Geod. 2000, 74, 128–133. [Google Scholar] [CrossRef]
  29. Groten, E. Parameters of Common Relevance of Astronomy, Geodesy, and Geodynamics. J. Geod. 2000, 74, 134–140. [Google Scholar] [CrossRef]
  30. Schönemann, P.H. A generalized solution of the orthogonal procrustes problem. Psychometrika 1966, 31, 1–10. [Google Scholar] [CrossRef]
  31. Watson, G.A. Computing Helmert transformations. J. Comput. Appl. Math. 2006, 197, 387–394. [Google Scholar] [CrossRef] [Green Version]
  32. Yang, J.; Wu, W.; Wu, Y.; Lian, J. Improved iterative calibration for triaxial accelerometers based on the optimal observation. Sensors 2012, 12, 8157–8175. [Google Scholar] [CrossRef] [Green Version]
  33. Panahandeh, G.; Skog, I.; Jansson, M. Calibration of the accelerometer triad of an inertial measurement unit, maximum likelihood estimation and Cramer-Rao bound. In Proceedings of the 2010 International Conference on Indoor Positioning and Indoor Navigation, Zurich, Switzerland, 15–17 September 2010; pp. 1–6. [Google Scholar] [CrossRef] [Green Version]
  34. Julier, S.J.; Uhlmann, J.K. New extension of the Kalman filter to nonlinear systems. In Proceedings of the AeroSense ’97, Orlando, FL, USA, 21–25 April 1997; Volume 3068. [Google Scholar] [CrossRef]
  35. Bar-Shalom, Y.; Li, X.; Kirubarajan, T. Estimation with Applications to Tracking and Navigation: Theory, Algorithms and Software; Wiley: Hoboken, NJ, USA, 2001. [Google Scholar]
  36. Dmitriev, V.I.; Ingtem, Z.G. Numerical differentiation using spline functions. Comput. Math. Model. 2012, 23, 312–318. [Google Scholar] [CrossRef]
  37. Roussel, C.; Verdun, J.; Cali, J.; Maia, M.; d’EU, J.F. Integration of a strapdown gravimeter system in an autonomous underwater vehicle. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2015, XL-5/W5, 199–206. [Google Scholar] [CrossRef] [Green Version]
  38. Maia, M. GRAVIMOB Cruise, L’Europe R/V; Sismer: Bologna, Italy, 2016. [Google Scholar] [CrossRef]
  39. Wessel, P.; Luis, J.F.; Uieda, L.; Scharroo, R.; Wobbe, F.; Smith, W.H.F.; Tian, D. The Generic Mapping Tools Version 6. Geochem. Geophys. Geosyst. 2019, 20, 5556–5564. [Google Scholar] [CrossRef] [Green Version]
  40. Ince, E.S.; Barthelmes, F.; Reißland, S.; Elger, K.; Förste, C.; Flechtner, F.; Schuh, H. ICGEM—15 years of successful collection and distribution of global gravitational models, associated services, and future plans. Earth Syst. Sci. Data 2019, 11, 647–674. [Google Scholar] [CrossRef] [Green Version]
  41. Lemoine, J.M.; Biancale, R.; Reinquin, F.; Bourgogne, S.; Gégout, P. CNES/GRGS RL04 Earth Gravity Field Models, from GRACE and SLR Data; GFZ Data Services: Potsdam, Germany, 2019. [Google Scholar] [CrossRef]
  42. Brockmann, J.M.; Schubert, T.; Schuh, W.D. An Improved Model of the Earth’s Static Gravity Field Solely Derived from Reprocessed GOCE Data. Surv. Geophys. 2021, 42, 277–316. [Google Scholar] [CrossRef]
  43. Zingerle, P.; Pail, R.; Gruber, T.; Oikonomidou, X. The combined global gravity field model XGM2019e. J. Geod. 2020, 94, 66. [Google Scholar] [CrossRef]
  44. Liang, W.; Li, J.; Xu, X.; Zhang, S.; Zhao, Y. A High-Resolution Earth’s Gravity Field Model SGG-UGM-2 from GOCE, GRACE, Satellite Altimetry, and EGM2008. Engineering 2020, 6, 860–878. [Google Scholar] [CrossRef]
  45. Verdun, J.; Damenet, N.; Cali, J. Moving-base vector gravimetry data processing based on optimal, physically sensible evolution models. In Proceedings of the IAG Symposium on Terrestrial Gravimetry: Static and Mobile Measurements (TG-SMM 2013), St. Petersburg, Russia, 1–4 October 2013; pp. 41–47. [Google Scholar]
  46. Newtonian Methods. In Numerical Optimization: Theoretical and Practical Aspects; Bonnans, J.F.; Gilbert, J.C.; Lemaréchal, C.; Sagastizábal, C.A. (Eds.) Springer: Berlin/Heidelberg, Germany, 2006; pp. 51–66. [Google Scholar] [CrossRef]
Figure 1. (a) QFLEX QA-3000 electrostatic accelerometers manufactured by Honeywell. (b) Picture showing GraviMob’s accelerometer sensor and its 24 bit digitiser. The two triads α and β each consisted of 3 accelerometers, the axes of which are mutually orthogonal and are shown in the red boxes. (c) Picture showing the glass sphere mounted in its frame itself screw-fixed into the AUV’s front case. (d) Overall view of AUV Aster x . The red arrow points to the front case where the sensor was installed.
Figure 1. (a) QFLEX QA-3000 electrostatic accelerometers manufactured by Honeywell. (b) Picture showing GraviMob’s accelerometer sensor and its 24 bit digitiser. The two triads α and β each consisted of 3 accelerometers, the axes of which are mutually orthogonal and are shown in the red boxes. (c) Picture showing the glass sphere mounted in its frame itself screw-fixed into the AUV’s front case. (d) Overall view of AUV Aster x . The red arrow points to the front case where the sensor was installed.
Remotesensing 14 02513 g001
Figure 2. View of IFREMER’s AUV Aster x . This submersible is 4.5 m long, 2.1 m wide (fins included) and 0.71 m in diameter with a total mass of 850 kg. It is powered by 50V/16kWh Li-ion batteries. Its autonomy is about 16 h. The maximum immersion depth of Aster x is 2650 m.
Figure 2. View of IFREMER’s AUV Aster x . This submersible is 4.5 m long, 2.1 m wide (fins included) and 0.71 m in diameter with a total mass of 850 kg. It is powered by 50V/16kWh Li-ion batteries. Its autonomy is about 16 h. The maximum immersion depth of Aster x is 2650 m.
Remotesensing 14 02513 g002
Figure 3. View of the AUV during the manipulation performed to adjust the internal parameters of the GraviMob system in one of the tilted positions. The AUV can thus be tilted along the roll, pitch and yaw axes.
Figure 3. View of the AUV during the manipulation performed to adjust the internal parameters of the GraviMob system in one of the tilted positions. The AUV can thus be tilted along the roll, pitch and yaw axes.
Remotesensing 14 02513 g003
Figure 4. View of the GraviMob system inside the climate chamber. The thumbnail shows an enlarged view of the additional electronic thermometer used for monitoring temperature variations inside the chamber.
Figure 4. View of the GraviMob system inside the climate chamber. The thumbnail shows an enlarged view of the additional electronic thermometer used for monitoring temperature variations inside the chamber.
Remotesensing 14 02513 g004
Figure 5. Graphs showing the variation of the scale factors k x (top), k y (middle) and k z (bottom) with temperature for the triad α . The red solid circles correspond to the values of the scale factors determined by our adjustment method and are plotted with their 1 sigma error bars. The blue curves represent the result of a polynomial fit with fourth-degree polynomials. Expressions of the polynomials are directly written out on the graphs. For all three polynomial fits, the coefficient of determination was found to be 0.987 ( k x ), 0.987 ( k y ) and 0.915 ( k z ). These coefficients did not show a significant increase for polynomial models of degrees greater than 4. The polynomial fits were validated by performing a chi-squared statistical test on the residuals with a 5% risk. This validation required the rejection of some observations that led to the failure of the chi-squared test. The same graphs are available for the triad β .
Figure 5. Graphs showing the variation of the scale factors k x (top), k y (middle) and k z (bottom) with temperature for the triad α . The red solid circles correspond to the values of the scale factors determined by our adjustment method and are plotted with their 1 sigma error bars. The blue curves represent the result of a polynomial fit with fourth-degree polynomials. Expressions of the polynomials are directly written out on the graphs. For all three polynomial fits, the coefficient of determination was found to be 0.987 ( k x ), 0.987 ( k y ) and 0.915 ( k z ). These coefficients did not show a significant increase for polynomial models of degrees greater than 4. The polynomial fits were validated by performing a chi-squared statistical test on the residuals with a 5% risk. This validation required the rejection of some observations that led to the failure of the chi-squared test. The same graphs are available for the triad β .
Remotesensing 14 02513 g005
Figure 6. Graphs showing the variation of the parameters τ x y (top left), τ y y (middle left), τ x z (top right), τ y z (middle right) and τ z z (bottom right) arranged as in the matrix (A15) with temperature for the triad α . The red solid circles correspond to the values of the parameters determined by our adjustment method and are plotted with their 1 sigma error bars. The blue curves represent the result of a polynomial fit with fourth-degree polynomials. Expressions of the polynomials are directly written out on the graphs. For all three polynomial fits, the coefficient of determination was found to be 0.988 ( τ x y ), 0.899 ( τ x z ), 0.975 ( τ y y ), 0.989 ( τ y z ) and 0.997 ( τ z z ). These coefficients did not show a significant increase for polynomial models of degrees greater than 4. The polynomial fits were validated by performing a chi-squared statistical test on the residuals with a 5% risk. This validation required the rejection of some observations that led to the failure of the chi-squared test. The same graphs are available for the triad β .
Figure 6. Graphs showing the variation of the parameters τ x y (top left), τ y y (middle left), τ x z (top right), τ y z (middle right) and τ z z (bottom right) arranged as in the matrix (A15) with temperature for the triad α . The red solid circles correspond to the values of the parameters determined by our adjustment method and are plotted with their 1 sigma error bars. The blue curves represent the result of a polynomial fit with fourth-degree polynomials. Expressions of the polynomials are directly written out on the graphs. For all three polynomial fits, the coefficient of determination was found to be 0.988 ( τ x y ), 0.899 ( τ x z ), 0.975 ( τ y y ), 0.989 ( τ y z ) and 0.997 ( τ z z ). These coefficients did not show a significant increase for polynomial models of degrees greater than 4. The polynomial fits were validated by performing a chi-squared statistical test on the residuals with a 5% risk. This validation required the rejection of some observations that led to the failure of the chi-squared test. The same graphs are available for the triad β .
Remotesensing 14 02513 g006
Figure 7. Graphs showing the variation of the accelerometer biases v x 0 (top), v y 0 (middle) and v z 0 (bottom) with temperature for the triad α . The red solid circles correspond to the values of the biases determined by our adjustment method and are plotted with their 1 sigma error bars. The blue curves represent the result of a polynomial fit with fourth-degree polynomials. Expressions of the polynomials are directly written out on the graphs. For all three polynomial fits, the coefficient of determination was found to be 0.993 ( v x 0 ), 0.996 ( v y 0 ) and 0.989 ( v z 0 ). These coefficients did not show a significant increase for polynomial models of degrees greater than 4. The polynomial fits were validated by performing a chi-squared statistical test on the residuals with a 5% risk. This validation required the rejection of some observations that led to the failure of the chi-squared test. The same graphs are available for the triad β .
Figure 7. Graphs showing the variation of the accelerometer biases v x 0 (top), v y 0 (middle) and v z 0 (bottom) with temperature for the triad α . The red solid circles correspond to the values of the biases determined by our adjustment method and are plotted with their 1 sigma error bars. The blue curves represent the result of a polynomial fit with fourth-degree polynomials. Expressions of the polynomials are directly written out on the graphs. For all three polynomial fits, the coefficient of determination was found to be 0.993 ( v x 0 ), 0.996 ( v y 0 ) and 0.989 ( v z 0 ). These coefficients did not show a significant increase for polynomial models of degrees greater than 4. The polynomial fits were validated by performing a chi-squared statistical test on the residuals with a 5% risk. This validation required the rejection of some observations that led to the failure of the chi-squared test. The same graphs are available for the triad β .
Remotesensing 14 02513 g007
Figure 8. Gravity survey in the Mediterranean Sea carried out with the GraviMob system and the AUV Aster x . The location of the eastern and western zones is indicated by the black boxes on map (b). Each zone comprises two profiles along which the surface gravity was measured by the SHOM. The underwater relief in these two zones is depicted on the bathymetric maps (a,c). Map (c) indicates the studied profile on which our data processing was tested. The profile was crossed twice following routes 8 and 9 respectively. Maps were made using Generic Mapping Tools software [39].
Figure 8. Gravity survey in the Mediterranean Sea carried out with the GraviMob system and the AUV Aster x . The location of the eastern and western zones is indicated by the black boxes on map (b). Each zone comprises two profiles along which the surface gravity was measured by the SHOM. The underwater relief in these two zones is depicted on the bathymetric maps (a,c). Map (c) indicates the studied profile on which our data processing was tested. The profile was crossed twice following routes 8 and 9 respectively. Maps were made using Generic Mapping Tools software [39].
Remotesensing 14 02513 g008
Figure 9. Temperature variation inside the sphere of the GraviMob system from the moment the AUV was launched until it was brought back on board. Once the temperature is stabilised, it remains remarkably constant at around 14.5 °C.
Figure 9. Temperature variation inside the sphere of the GraviMob system from the moment the AUV was launched until it was brought back on board. Once the temperature is stabilised, it remains remarkably constant at around 14.5 °C.
Remotesensing 14 02513 g009
Figure 10. Depth of the underwater relief and surface free air gravity anomaly both provided by the SHOM as a function of distance travelled along the studied profile.
Figure 10. Depth of the underwater relief and surface free air gravity anomaly both provided by the SHOM as a function of distance travelled along the studied profile.
Remotesensing 14 02513 g010
Figure 11. Variations of the gravity vector components g e , g n , g u along the studied profile along route 8 (solid blue curve) and route 9 (solid green curve), calculated according to our UKF-based processing method. The results of route 9 have been reversed to be comparable with the results of route 8.
Figure 11. Variations of the gravity vector components g e , g n , g u along the studied profile along route 8 (solid blue curve) and route 9 (solid green curve), calculated according to our UKF-based processing method. The results of route 9 have been reversed to be comparable with the results of route 8.
Remotesensing 14 02513 g011
Figure 12. Variations of the formal errors on the east (blue solid curve), north (green solid curve) and vertical (red solid curve) gravity field components obtained from the covariance matrix given by the unscented Kalman filter. The three errors all stabilise after about 1000 time steps, which corresponds to a distance travelled on the profile of 750 m. Once stabilised, the value of σ g u is almost four times smaller than that of σ g e and σ g n .
Figure 12. Variations of the formal errors on the east (blue solid curve), north (green solid curve) and vertical (red solid curve) gravity field components obtained from the covariance matrix given by the unscented Kalman filter. The three errors all stabilise after about 1000 time steps, which corresponds to a distance travelled on the profile of 750 m. Once stabilised, the value of σ g u is almost four times smaller than that of σ g e and σ g n .
Remotesensing 14 02513 g012
Figure 13. Variations of gravity anomaly δ g on routes 8 (blue solid curve) and 9 (green solid curve) respectively, calculated from the gravity vector intensity g ¯ n once high-frequency variations were smoothed out using a 3000 m moving window filter.
Figure 13. Variations of gravity anomaly δ g on routes 8 (blue solid curve) and 9 (green solid curve) respectively, calculated from the gravity vector intensity g ¯ n once high-frequency variations were smoothed out using a 3000 m moving window filter.
Remotesensing 14 02513 g013
Figure 14. Variations of gravity anomaly trends deduced from SHOM (red solid curve) and GraviMob measurements on routes 8 (blue solid curve) and 9 (green solid curve). These trends were obtained by fitting second degree polynomials.
Figure 14. Variations of gravity anomaly trends deduced from SHOM (red solid curve) and GraviMob measurements on routes 8 (blue solid curve) and 9 (green solid curve). These trends were obtained by fitting second degree polynomials.
Remotesensing 14 02513 g014
Figure 15. Left-hand plot: variations of the surface gravity anomaly deduced from that measured on route 9 ( δ g 9 ) by upward continuation using transfer function f 8 . The SHOM anomaly is represented by the red solid curve. Right-hand plot: histogram of the differences between the two anomalies at each measured point of the profile.
Figure 15. Left-hand plot: variations of the surface gravity anomaly deduced from that measured on route 9 ( δ g 9 ) by upward continuation using transfer function f 8 . The SHOM anomaly is represented by the red solid curve. Right-hand plot: histogram of the differences between the two anomalies at each measured point of the profile.
Remotesensing 14 02513 g015
Figure 16. Left-hand plot: variations of the surface gravity anomaly deduced from that measured on route 8 ( δ g 8 ) by upward continuation using transfer function f 9 . The SHOM anomaly is represented by the red solid curve. Right-hand plot: histogram of the differences between the two anomalies at each measured point of the profile.
Figure 16. Left-hand plot: variations of the surface gravity anomaly deduced from that measured on route 8 ( δ g 8 ) by upward continuation using transfer function f 9 . The SHOM anomaly is represented by the red solid curve. Right-hand plot: histogram of the differences between the two anomalies at each measured point of the profile.
Remotesensing 14 02513 g016
Figure 17. (a) Parallelepiped triad versus (b) pyramid-shaped triad. It should be noted that the pyramid-shaped triad allows equal wear of each accelerometer.
Figure 17. (a) Parallelepiped triad versus (b) pyramid-shaped triad. It should be noted that the pyramid-shaped triad allows equal wear of each accelerometer.
Remotesensing 14 02513 g017
Table 1. Technical data of QFLEX QA-3000 electrostatic accelerometers. The letter “ g ” stands for an acceleration of 10 m s 2 , that is 10 6 mGal. The possible measuring range of the accelerometers has been deliberately narrowed so as to lessen the digitizing error. A voltage offset has therefore to be applied to the accelerometer of each triad, the sensitive axis of which is most often in a direction close to the vertical. The electrical arrangement of the sensor only allows for specific force intensities below g/5.32 to be measured. Above this value, the accelerometers enter saturation plage. Therefore, the value of the limit tilt angle beyond which horizontal accelerometers enter saturation is 10.8 deg, value deduced from arcsin ( 1 / 5.32 ) . A voltage source that offers at least 5 μ V accuracy is crucial since a variation of such a value would be interpreted as a 1 mGal variation of the acceleration.
Table 1. Technical data of QFLEX QA-3000 electrostatic accelerometers. The letter “ g ” stands for an acceleration of 10 m s 2 , that is 10 6 mGal. The possible measuring range of the accelerometers has been deliberately narrowed so as to lessen the digitizing error. A voltage offset has therefore to be applied to the accelerometer of each triad, the sensitive axis of which is most often in a direction close to the vertical. The electrical arrangement of the sensor only allows for specific force intensities below g/5.32 to be measured. Above this value, the accelerometers enter saturation plage. Therefore, the value of the limit tilt angle beyond which horizontal accelerometers enter saturation is 10.8 deg, value deduced from arcsin ( 1 / 5.32 ) . A voltage source that offers at least 5 μ V accuracy is crucial since a variation of such a value would be interpreted as a 1 mGal variation of the acceleration.
Technical FeatureValue or Range of Values
Possible measuring range ± 60 g ( ± 60 × 10 6   mGal )
Tailored measuring range ± 1 5.32 g ( ± 187 , 970 mGal )
Bias< 4 × 10 3 g ( 4000 mGal )
Scale factor (output current)between 1.20 and 1.46   mA / g
Scale factor (output voltage)between 4.80 and 5.84   V / g
Temperature sensitivity15 μ g / C (15 mGal / C )
Table 2. Mean and standard deviation of the residuals (Equation (35)) obtained with our adjustment method and corresponding values of the internal parameters (scale factors, misalignment components, biases) for the triads α and β .
Table 2. Mean and standard deviation of the residuals (Equation (35)) obtained with our adjustment method and corresponding values of the internal parameters (scale factors, misalignment components, biases) for the triads α and β .
ResidualsTriad α Triad β Unit
Mean00mGal
Standard 1.4 3.1 mGal
ParametersTriad α Triad β Unit
(value ± 1 σ error)
k x 5.3903 ± 3 × 10 4 5.3484 ± 2 × 10 4 μ V mGal 1
k y 5.3892 ± 3 × 10 4 5.3426 ± 2 × 10 4 μ V mGal 1
k z 5.489 ± 1 × 10 3 5.2792 ± 6 × 10 4 μ V mGal 1
τ x y + 0.0003216 ± 2 × 10 7 + 0.0008731 ± 2 × 10 7
τ x z + 0.003873 ± 4 × 10 6 0.000903 ± 3 × 10 6
τ y y + 1.0000000000 ± 1 × 10 10 + 1.0000003811 ± 2 × 10 10
τ y z 0.002580 ± 4 × 10 6 0.002332 ± 3 × 10 6
τ z z + 1.00001083 ± 3 × 10 8 + 1.000003126 ± 3 × 10 9
v x 0 14.89 ± 0.02 1.62 ± 0.02 mV
v y 0 9.12610 ± 2 × 10 5 14.56 ± 0.02 mV
v z 0 73.3 ± 0.8 31.3 ± 0.4 mV
Table 3. Euler angle values that defined the matrices C s o b for triads α and β respectively. The specific force measurements needed to determine these values were acquired for N = 31 tilted positions of the AUV. During this experiment, the average temperature inside the laboratory was 12.4 °C, and the temperature variation between the beginning and end of the experiment was less than 2 °C. The differences between the orientations of triads α and β are mainly due to the mounting of the sphere inside the AUV. It is clear that the values of these angles must be re-estimated regularly and at each new installation of the sphere in the AUV. These results indicate that the axes of triad α are almost parallel to those of triad β , which is consistent with the assembly of the two triads shown in Figure 1b.
Table 3. Euler angle values that defined the matrices C s o b for triads α and β respectively. The specific force measurements needed to determine these values were acquired for N = 31 tilted positions of the AUV. During this experiment, the average temperature inside the laboratory was 12.4 °C, and the temperature variation between the beginning and end of the experiment was less than 2 °C. The differences between the orientations of triads α and β are mainly due to the mounting of the sphere inside the AUV. It is clear that the values of these angles must be re-estimated regularly and at each new installation of the sphere in the AUV. These results indicate that the axes of triad α are almost parallel to those of triad β , which is consistent with the assembly of the two triads shown in Figure 1b.
AnglesTriad α Triad β Unit
θ x 0.310 0.297 decimal degrees
θ y 0.153 0.192 decimal degrees
θ z 177.334 177.142 decimal degrees
Table 4. Values of the standard deviations that enter the expression of the covariance matrix Q of the evolution model (49).
Table 4. Values of the standard deviations that enter the expression of the covariance matrix Q of the evolution model (49).
PositionValueUnit
σ λ ¨ 5 × 10 6 deg s 2
σ φ ¨ 4 × 10 6 deg s 2
σ h ¨ 0.1 m s 2
AttitudeValueUnit
σ δ ¨ 0.8 deg s 2
σ χ ¨ 0.5 deg s 2
σ η ¨ 1.7 deg s 2
GravityValueUnit
σ g ¨ e 1 × 10 3 mGal s 2
σ g ¨ n 1 × 10 3 mGal s 2
σ g ¨ u 1 × 10 3 mGal s 2
Table 5. Properties of recent global gravity field models used to calculate estimates of the vertical deflection components on the study profile. The column entitled Data indicates by specific letters which datasets have been used to compute the geopotential model in question, i.e., A for satellite altimetry data, S for satellite gravity data (GRACE, GOCE, LAGEOS), G for ground data (terrestrial, shipborne and airborne gravity measurements) and T for topography data. The letter M represents a model that derives from a previous model taken as an a priori. The first two models (no 1 and no 2) calculated from satellite data only have a spatial resolution of 67 km. The last two (no 3 and no 4) assimilate a wide variety of gravity data, which allows for a finer spatial resolution of 9 km.
Table 5. Properties of recent global gravity field models used to calculate estimates of the vertical deflection components on the study profile. The column entitled Data indicates by specific letters which datasets have been used to compute the geopotential model in question, i.e., A for satellite altimetry data, S for satellite gravity data (GRACE, GOCE, LAGEOS), G for ground data (terrestrial, shipborne and airborne gravity measurements) and T for topography data. The letter M represents a model that derives from a previous model taken as an a priori. The first two models (no 1 and no 2) calculated from satellite data only have a spatial resolution of 67 km. The last two (no 3 and no 4) assimilate a wide variety of gravity data, which allows for a finer spatial resolution of 9 km.
Model NumberModel NameYearDataDegreeReference
1EIGEN-GRGS.RL04.MEAN-FIELD2019S300[41]
2GO_CONS_GCF_2_TIM_R62019S300[42]
3XGM2019e_21592019A, G, S, T2190[43]
4SGG-UGM-22020A, M, S2190[44]
Table 6. Values of statistical indicators (minimum, maximum, mean and standard deviation) on the components of the vertical deviation calculated from 1000 measurement points on the studied profile. The number of significant figures was set according to the value of the standard deviation. The 4 models were found to be consistent with each other, with the exception of the north–south deflection components calculated from models 1 and 2. The values of the vertical deviation were calculated at the height of the profile studied with respect to the reference ellipsoid estimated at -1850 m. This estimate was determined by considering a depth of 1900 m below the geoid surface, which in this area is on average 50 m above the surface of the reference ellipsoid. The mean values derived from these models represent reasonable estimates of the vertical deflection components on the profile studied.
Table 6. Values of statistical indicators (minimum, maximum, mean and standard deviation) on the components of the vertical deviation calculated from 1000 measurement points on the studied profile. The number of significant figures was set according to the value of the standard deviation. The 4 models were found to be consistent with each other, with the exception of the north–south deflection components calculated from models 1 and 2. The values of the vertical deviation were calculated at the height of the profile studied with respect to the reference ellipsoid estimated at -1850 m. This estimate was determined by considering a depth of 1900 m below the geoid surface, which in this area is on average 50 m above the surface of the reference ellipsoid. The mean values derived from these models represent reasonable estimates of the vertical deflection components on the profile studied.
EW Vertical Deflection η (Arcsec)NS Vertical Deflection ξ (Arcsec)
Model1234Model1234
Mean3.563.533.63.0Mean—8.43—8.73—9—10
Min3.463.302.92.7Min—8.46—8.77—13—15
Max3.643.734.74.0Max—8.35—8.64—5—6
Std0.050.120.50.4Std0.030.0433
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Verdun, J.; Roussel, C.; Cali, J.; Maia, M.; D’Eu, J.-F.; Kharbou, O.; Poitou, C.; Ammann, J.; Durand, F.; Bouhier, M.-É. Development of a Lightweight Inertial Gravimeter for Use on Board an Autonomous Underwater Vehicle: Measurement Principle, System Design and Sea Trial Mission. Remote Sens. 2022, 14, 2513. https://doi.org/10.3390/rs14112513

AMA Style

Verdun J, Roussel C, Cali J, Maia M, D’Eu J-F, Kharbou O, Poitou C, Ammann J, Durand F, Bouhier M-É. Development of a Lightweight Inertial Gravimeter for Use on Board an Autonomous Underwater Vehicle: Measurement Principle, System Design and Sea Trial Mission. Remote Sensing. 2022; 14(11):2513. https://doi.org/10.3390/rs14112513

Chicago/Turabian Style

Verdun, Jérôme, Clément Roussel, José Cali, Marcia Maia, Jean-François D’Eu, Ossama Kharbou, Charles Poitou, Jérôme Ammann, Frédéric Durand, and Marie-Édith Bouhier. 2022. "Development of a Lightweight Inertial Gravimeter for Use on Board an Autonomous Underwater Vehicle: Measurement Principle, System Design and Sea Trial Mission" Remote Sensing 14, no. 11: 2513. https://doi.org/10.3390/rs14112513

APA Style

Verdun, J., Roussel, C., Cali, J., Maia, M., D’Eu, J. -F., Kharbou, O., Poitou, C., Ammann, J., Durand, F., & Bouhier, M. -É. (2022). Development of a Lightweight Inertial Gravimeter for Use on Board an Autonomous Underwater Vehicle: Measurement Principle, System Design and Sea Trial Mission. Remote Sensing, 14(11), 2513. https://doi.org/10.3390/rs14112513

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop