2.1. System Formulation

The DVL uses the principle that by reflecting three or more non-coplanar radio/sound beams off a surface and measuring the Doppler shifts, the velocity of the body with respect to that surface can be obtained (Groves, Reference Groves2008). Most systems use the Janus configuration with four beams. The surface-referenced velocity can be obtained as

(1)$${\bf v}^d = {\bf K\Delta f}$$

The superscript d denotes the DVL's coordinate frame D naturally defined by the beam spatial configuration. The vector Δf is formed by Doppler shifts along each beam and the matrix K depends on the transmitted sound wave frequency, the spatial configuration of beams and the sound speed in water. The first two factors of K are fixed for a DVL unit, but the sound speed may vary with temperature, depth and salinity by a few percent (Groves, Reference Groves2008). We rewrite the above relationship as

(2)$${\bf v}^d = {\bf K\Delta f} = \displaystyle{1 \over k}{\bf K}_s {\bf \Delta f} = \displaystyle{1 \over k}{\bf y}_{DVL} \Leftrightarrow {\bf y}_{DVL} = k{\bf v}^d $$

where y_DVL ≜ K_sΔf denotes the DVL output, k is a scalar factor that accounts for the change of sound velocity, and K_s is a scaled version of K that is nearly invariant to water conditions.

In contrast to the DVL directly providing the velocity relative to the seabed (bottom-lock), the IMU measures the angular velocity and non-gravitational acceleration with respect to the inertial space. Numerical integrations must be carried out to derive attitude, velocity or position information, which is known as the inertial navigation computation procedure (Groves, Reference Groves2008; Wu and Pan, Reference Wu and Pan2013b). As gyroscopes and accelerometers are subject to errors like bias and noise, the computed inertial navigation result is prone to error drift accumulating with time.

Denote by N the local level reference frame, by B the IMU body frame, by I the inertial non-rotating frame, and by E the Earth frame. The navigation (attitude, velocity and position) rate equations in the reference n-frame are well known as (Titterton and Weston, Reference Titterton and Weston2004; Groves, Reference Groves2008; Savage, Reference Savage2007)

(3)$$\dot{\bf C}_b^n = {\bf C}_b^n \left( {{\bf \omega} _{nb}^b \times} \right){\comma}\quad {\bf \omega} _{nb}^b = {\bf \omega} _{ib}^b - {\bf b}_g - {\bf C}_n^b \left( {{\bf \omega} _{ie}^n + {\bf \omega} _{en}^n} \right)\comma$$

(4)$$\dot{\bf v}^n = {\bf C}_b^n \left( {{\bf f}^b - {\bf b}_a} \right) - \left( {2{\bf \omega} _{ie}^n + {\bf \omega} _{en}^n} \right) \times {\bf v}^n + {\bf g}^n $$

(5)$$\dot{\bf p} = {\bf R}_c {\bf v}^n $$

where ${\bf C}_b^n $ denotes the attitude matrix from the body frame to the reference frame, vⁿ the velocity relative to the Earth, ${\bf \omega} _{ib}^b $ the error-contaminated body angular rate measured by gyroscopes in the body frame, f^b the error-contaminated specific force measured by accelerometers in the body frame, ${\bf \omega} _{ie}^n $ the Earth rotation rate with respect to the inertial frame, ${\bf \omega} _{en}^n $ the angular rate of the reference frame with respect to the Earth frame, ${\bf \omega} _{nb}^b $ the body angular rate with respect to the reference frame, and gⁿ the gravity vector. The skew symmetric matrix (·×) is defined so that the cross product satisfies a × b = (a×)b for two arbitrary vectors. The gyroscope bias b_g and the accelerometer bias b_a are taken into consideration as approximately random constants, i.e., ${\dot{\bf b}_g} = {\dot{\bf b}_a} = {\bf 0}$. The position ${\bf p} \mathop{=}\limits^\Delta \left[ {\matrix{\lambda & L & h \cr}} \right]^T $ is described by the angular orientation of the reference frame relative to the Earth frame, commonly expressed as longitude λ, latitude L and height h above the Earth's surface. In the context of a specific local level frame choice, e.g., North-Up-East, ${\bf v}^n \mathop{=}\limits^\Delta \left[ {\matrix{{v_N} & {v_U} & {v_E}\cr}} \right]^T $ and the local curvature matrix is explicitly expressed as

(6)$${\bf R}_c = \left[ {\matrix{ 0 & 0 & {\displaystyle{1 \over {(R_E + h)\cos L}}} \cr {\displaystyle{1 \over {R_N + h}}} & 0 & 0 \cr 0 & 1 & 0 \cr}} \right]$$

where R _E and R _N are respectively the transverse radius of curvature and the meridian radius of curvature of the reference ellipsoid.

Using Equation (2), the derived velocity from IMU is related to the DVL output by

(7)$${\bf y}_{DVL} = k{\bf C}_b^d {\bf C}_n^b {\bf v}^n $$

where ${\bf C}_b^d $ is the misalignment attitude matrix of the d-frame with respect to the b-frame. This misalignment matrix and the scale factor are both regarded as random constants. As with previous studies (Hegrenæs and Berglund, Reference Hegren, S. and Berglund2009; Joyce, Reference Joyce1989; Kinsey and Whitcomb, Reference Kinsey and Whitcomb2006; Reference Kinsey and Whitcomb2007; Tang et al., Reference Tang, Wang, Li and Wu2013; Troni et al., Reference Troni, Kinsey, Yoerger and Whitcomb2012; Troni and Whitcomb, Reference Troni and Whitcomb2010), the translational misalignment between DVL and IMU will not be considered hereafter.

We see from Equations (3), (4) and (7) that in addition to attitude/velocity/position that are of immediate interest to us, the combined IMU/DVL navigation also necessitates the finding of such parameters as inertial sensor biases (b_g and b_a), DVL scale factor k and misalignment matrix ${\bf C}_b^d $. Insufficient knowledge of these parameters could result in degrading the combined navigation accuracy.

2.2. State Observability Analysis

Without any other external sensors, is it possible to determine the above parameters? From a viewpoint of a control system, this question relates to the (global) observability of the system state (Chen, Reference Chen1999; Wu et al., Reference Wu, Zhang, Wu, Hu and Hu2012). In such a case, the system model is given by Equations (3)–(5) and the observation model is given by Equation (7), with the IMU measurement as the system input and the DVL measurement as the system output.

Hereafter we use y to replace y_DVL for notational brevity. From Equation (7), we have ${\bf v}^n = {\bf C}_b^n {\bf C}_d^b {\bf y}/k$. Substituting into Equation (4) and using Equation (3) yield

(8)$${{\left( {{\bf C}_b^n \left( {{\bf \omega}_{nb}^b \times} \right){\bf C}_d^b {\bf y} + {\bf C}_b^n {\bf C}_d^b \dot{\bf y}} \right)} / k} = {\bf C}_b^n \left( {{\bf f}^b - {\bf b}_a} \right) - \left( {2{\bf \omega}_{ie}^n + {\bf \omega}_{en}^n} \right) \times {{\left( {{\bf C}_b^n {\bf C}_d^b {\bf y}} \right)} / k} + {\bf g}^n $$

or equivalently,

(9)$$k\left( {{\bf C}_n^b {\bf g}^n + {\bf f}^b - {\bf b}_a} \right) = \left( {\left( {{\bf \omega} _{ie}^b + {\bf \omega}_{ib}^b - {\bf b}_g} \right) \times} \right){\bf C}_d^b {\bf y} + {\bf C}_d^b \dot{\bf y}$$

Troni et al. (Reference Troni, Kinsey, Yoerger and Whitcomb2012) and Troni and Whitcomb (Reference Troni and Whitcomb2010) used a quite coarse simplification of this equation for DVL calibration, for instance, neglecting the term ${\bf C}_n^b {\bf g}^n $.

On any trajectory segment with constant ${\bf C}_n^b $, the quantities ${\bf \omega} _{ie}^b $ and ${\bf \omega} _{ib}^b - {\bf b}_g $ keep almost unchanged with the moderate speed of an underwater vehicle. As the vector magnitude $\Vert {{\bf \omega} _{ib}^b - {\bf b}_g} \Vert \approx \Vert {{\bf \omega} _{ie}^b} \Vert \approx 7\hskip-2pt \cdot \hskip-2pt 3 \times 10^{ - 5} \,{\rm rad}/{\rm s}$ is very small, Equation (9) is reasonably approximated by

(10)$$k\left( {{\bf C}_n^b {\bf g}^n + {\bf f}^b - {\bf b}_a} \right) = {\bf C}_d^b \dot{\bf y}$$

Time derivative of the above equality is

(11)$$k\dot{\bf f}^b = {\bf C}_d^b \ddot{\bf y}$$

from which we obtain $k = \pm {{\Vert {\ddot{\bf y}} \Vert} / {\Vert {{{\dot{\bf f}}^b}} \Vert}}$ whenever $\Vert {{{\dot{\bf f}}^b}} \Vert$ is non-vanishing, as an attitude matrix does not change the magnitude of a vector. It is trivial to solve the sign ambiguity. For example, if the b-frame is roughly aligned with the d-frame, which is often the case in practice, the positive sign will obviously be the right option. So the DVL scale factor is now a known quantity. Note that Equation (11) is valid for any segment of this type, so if there exist two such segments that ${\dot{\bf f}^b}$ (or equivalently $$\ddot{\bf y}$$) have different directions, the misalignment matrix ${\bf C}_d^b $ can be determined according to Lemma 1 in the Appendix.

On any trajectory segment with fast-changing ${\bf C}_n^b $, the quantity ${\bf \omega} _{ib}^b $ is much larger in magnitude than ${\bf \omega} _{ie}^b $ or b_g (as far as a quality IMU is concerned). Equation (9) can be approximated as

(12)$${\bf C}_n^b {\bf g}^n = {{\left( {\left( {{\bf \omega}_{ib}^b \times} \right){\bf C}_d^b {\bf y} + {\bf C}_d^b \dot{\bf y}} \right)} / {k - {\bf f}^b + {\bf b}_a}} $$

Taking the norm on both sides gives

(13)$$g = \Vert {{\bf g}^n} \Vert = \Vert {{\bf \alpha} + {\bf b}_a} \Vert $$

which is a quadratic equation on the accelerometer bias b_a with ${\bf \alpha} \mathop{=}\limits^\Delta {{\left( {\left( {{\bf \omega} _{ib}^b \times} \right){\bf C}_d^b {\bf y}\; + {\bf C}_d^b \dot{\bf y}} \right)} \big/ {k - {\bf f}^b}} $. According to Lemma 2 in the Appendix, we know that b_a will be determined if the vectors α, at all times on segments of this type, are non-coplanar, or ∑αα^T is non-singular (Wu et al., Reference Wu, Zhang, Wu, Hu and Hu2012). This requirement is naturally met for practical turning in water, as shown in the simulation section.

Then by the chain rule of the attitude matrix, ${\bf C}_n^b $ at any time on this segment satisfies

(14)$${\bf C}_n^b = {\bf C}_{b(t_0 )}^{b(t)} {\bf C}_{n(t_0 )}^{b(t_0 )} {\bf C}_{n(t)}^{n(t_0 )} \mathop{=}\limits^\Delta {\bf C}_{b(t_0 )}^{b(t)} {\bf C}_n^b (t_0 ){\bf C}_{n(t)}^{n(t_0 )} $$

where ${\bf C}_n^b (t_0 )$ denotes the initial attitude matrix at the beginning of this segment, and ${\bf C}_{b(t_0 )}^{b(t)} $ and ${\bf C}_{n(t_0 )}^{n(t)} $, respectively, encode the attitude changes of the body frame and the reference frame. Substituting Equation (14) into Equation (12),

(15)$${\bf C}_n^b (t_0 ){\bf C}_{n(t)}^{n(t_0 )} {\bf g}^n = {\bf C}_{b(t)}^{b(t_0 )} ({\bf \alpha} + {\bf b}_a )$$

${\bf C}_n^b (t_0 )$ is solvable as there always exist two time instants that ${\bf C}_{n(t)}^{n(t_0 )} {\bf g}^n $ have different directions due to the Earth rotation (Wu et al., Reference Wu, Zhang, Wu, Hu and Hu2012). Also solvable is the attitude matrix ${\bf C}_n^b $ using Equation (14). Then the velocity vⁿ can be determined by Equation (7), and the gyroscope bias b_g can be determined by Equation (3).

The above analysis is summarised in the theorem below.

Theorem 1 (State Observability): If the vehicle trajectory contains segments of constant attitude with linearly independent ${\dot{\bf f}^b}$ (Type-I), as well as turning segments on which the vectors α as defined in Equation (13) are non-coplanar (Type-II), then the system of Equations (3)–(5) and (7) is observable in attitude, velocity, inertial sensor biases, and DVL scale factor and misalignment matrix.

Here are some explanations on the condition of linearly independent ${\dot{\bf f}^b}$, the rate of the specific force (sum of external forces except gravitation).

Remark 1: For a multiple-thruster-propelled underwater vehicle as in Kinsey and Whitcomb (Reference Kinsey and Whitcomb2006; Reference Kinsey and Whitcomb2007), it is not difficult to fulfil this condition, for example by thruster switching. This will normally create driving forces in different directions in the IMU body frame, resulting in linearly independent ${\dot{\bf f}^b}$.

Remark 2: For an underwater vehicle with a single thruster, it is tricky to fulfil this condition while keeping constant attitude. No matter the thruster force or the water resistance, its direction is fixed relative to IMU or DVL. If the vehicle moves at the same depth, linearly independent ${\dot{\bf f}^b}$ might only occur instantaneously at both start and end times of Type-I segments (a situation much like that in Tang et al. (Reference Tang, Wu, Wu, Wu, Hu and Shen2009)). This excitation is in practice not sufficient to produce a good estimate. Consider an example where $\dot{\bf f}^b = \left[ {\matrix{{\dot f_x^b} & 0 & 0 \cr}} \right]^T $ for all Type-I segments, which means the thrusters' force aligns with x-axis of the IMU body frame. Suppose the rotation sequence from d-frame to b-frame is first around y-axis (yaw angle, ψ), followed by z-axis (pitch angle, θ) and then by x-axis (roll angle, φ), ${\bf C}_d^b $ can be re-parameterised in Euler angles as

(16)$${\bf C}_d^b = \left[ {\matrix{ {\cos \theta \cos \psi} & {\sin \theta} & { - \cos \theta \sin \psi} \cr {\sin \phi \sin \psi - \cos \phi \cos \psi \sin \theta} & {\cos \phi \cos \theta} & {\cos \psi \sin \phi + \cos \phi \sin \theta \sin \psi} \cr {\cos \phi \sin \psi + \cos \psi \sin \phi \sin \theta} & { - \cos \theta \sin \phi} & {\cos \phi \cos \psi - \sin \phi \sin \theta \sin \psi} \cr}} \right].$$

Substituting into Equation (11)

(17)$$k\,{\dot f_x^b} \left[ {\matrix{ {\cos \theta \cos \psi} & {\sin \theta} & { - \cos \theta \sin \psi} \cr}} \right]^T = \ddot{\bf y}$$

which is irrelevant to the roll angle φ around x-axis. The other two angles can be computed by element comparison, while the angle φ is not estimable. This case is right what we encountered in land vehicle navigation subject to the non-holonomic constraint (Wu et al., Reference Wu, Wu, Hu and Hu2009). It is the angle around the thruster force that is inestimable, so if the thruster force was in a general direction, it can be imagined that all three Euler angles would be affected.

Remark 3: Even if ${\dot{\bf f}^b}$ (or equivalently $\ddot{\bf y}$) on the segments of constant attitude are linearly dependent, all states in Theorem 1 are still observable except the angle around the direction of ${\dot{\bf f}^b}$. This interesting conclusion can be obtained from the proof of Theorem 1. We only need to show that the products ${\bf C}_d^b {\bf y}$ and ${\bf C}_d^b \dot{\bf y}$ are known quantities in this scenario, in spite of the undetermined ${\bf C}_d^b $. Denote by η the unit direction of linearly dependent $\ddot{\bf y}$. For a single-thruster underwater vehicle, both the DVL velocity y and the velocity rate $\dot{\bf y}$ almost certainly align with the direction of η. From Equation (11), for some time t ₁ on this segment $k\dot{\bf f}^b (t_1 ) = {\bf C}_d^b \ddot{\bf y}(t_1 ) = \Vert {\ddot{\bf y}(t_1 )} \Vert{\bf C}_d^b {\bf \eta} $ where $\ddot{\bf y}({t_1})$ is non-zero. Therefore, ${\bf C}_d^b {\bf y} = {{\Vert {\bf y} \Vert{\bf C}_d^b {\bf \eta} = k\Vert {\bf y} \Vert\dot{\bf f}^b (t_1 )} / {\Vert {\ddot{\bf y}(t_1 )} \Vert}}$ and ${\bf C}_d^b \dot{\bf y} = {{\Vert {\dot{\bf y}} \Vert{\bf C}_d^b {\bf \eta} = k\Vert {\dot{\bf y}} \Vert\dot{\bf f}^b (t_1 )} / {\Vert {\ddot{\bf y}(t_1 )} \Vert}}$ are both functions of known quantities.

Remark 4: Ascending or descending motion makes it possible for a single-thruster underwater vehicle to have linearly independent ${\dot{\bf f}^b}$. The water resistance due to ascending/descending velocity change will incur non-zero ${\dot{\bf f}^b}$ that deviates in direction from the thruster force. The ascending/descending motions are realised by unloading/loading the ballast tank.

2.3. Derived Ideal Observers

The above observability analysis not only gives a ‘yes or no’ answer to state estimability, but provides us valuable insights on designing observers. This is an additional advantage of global observability analysis (Wu et al., Reference Wu, Zhang, Wu, Hu and Hu2012). Specifically, Equation (11) allows us to construct an ideal observer to estimate DVL parameters, using IMU and DVL measurements on Type-I segments. No additional external sensors like GPS or LBL beacons are required. Suppose the time interval $\left[ {\matrix{{t_s^j} & {t_e^j} \cr}} \right]$ corresponds to the j-th Type-I segment (j = 1,2,…). For any $t \in \left[ {\matrix{{t_s^j} & {t_e^j}\cr}} \right]$, integrating Equation (11) twice over the subinterval $\left[ {\matrix{{t_s^j} & t \cr}} \right]$

(18)$$k{\bf \beta} (t) = {\bf C}_d^b {\bf \gamma} (t)$$

where ${\bf \beta} (t) \mathop{=}\limits^\Delta \int_{t_s^j} ^t {{\bf f}^b} d\tau - {\bf f}^b \left( {t_s^j} \right)\left( {t - t_s^j} \right)$ and ${\bf \gamma} (t) \mathop{=}\limits^\Delta {\bf y}\left( t \right) - {\bf y}\left( {t_s^j} \right) - \dot{\bf y}\left( {t_s^j} \right)\left( {t - t_s^j} \right)$. Compared with Equation (11), this integral form spares the differentiation of the measurements that would be subject to noise amplification (Wu and Pan, Reference Wu and Pan2013a; Wu et al., Reference Wu, Wang and Hu2014). The above equation obviously applies to all Type-I segments. Denote $\Omega = \bigcup\limits_j {\left[ {\matrix{{t_s^{\;j}} & {t_e^{\hskip2pt j}} \cr}} \right]} $, we have the following observer to calibrate the DVL parameters.

Ideal Observer – DVL Calibration (IO-DVLC):

1. 1) Compute $k = {{\Vert {{\bf \gamma }(t)} \Vert} / {\Vert {{\bf \beta }(t)} \Vert}}$ for each t ∈ Ω;

2. 2) Obtain ${\bf C}_d^b $ by solving the constrained optimization problem $\mathop {\min}\limits_{{\bf C} \in SO(3)} \int_\Omega {{{\Vert {k{\bf \beta }(t) - {\bf C\gamma }(t)} \Vert}^2}dt} $ as done in (Wu and Pan, Reference Wu and Pan2013a).

The “ideal observer” will be used to verify the analytic conclusions above.