4.1.1. Recovery of angular orientation

The tracking algorithm essentially relies on the knowledge of the angular orientation of the sensor system at every time step. The angular orientation at time t is uniquely determined by an orthogonal matrix

(1)$${\bf R}( t) = \left[\matrix{{\bf e}^1_x( t) & {\bf e}^1_y( t) & {\bf e}^1_z( t) }\right]$$

establishing an active rotation of the inertial (global) coordinate axes $( {\bf e}^0_x,\; \, {\bf e}^0_y,\; \, {\bf e}^0_z)$ onto the sensor-fixed (local) coordinate frame (¹e_x,  ¹e_y,  ¹e_z), see Figure 1. The kinematic differential equation governing the evolution of R(t) reads (Wittenburg, Reference Wittenburg2008)

(2)$$\dot{\bf R}( t) = {}^0{\bf \Omega}( t) \, {\bf R}( t) ,\; \quad {\bf R}( 0) = {\bf R}_0$$

with

(3)$${}^0{\bf \Omega} = \left[\matrix{0 & -{}^0\omega_z & {}^0\omega_y \cr {}^0\omega_z & 0 & -{}^0\omega_x \cr - {}^0\omega_y & {}^0\omega_x & 0 }\right]$$

being the skew-symmetric 3 × 3 matrix corresponding to the vector ⁰ω = [ ⁰ω _x ⁰ω _y ⁰ω _z ]^T of angular velocity components. To distinguish between vector components related to the corotating sensor frame and those related to the inertial frame, a left-hand superscript index is introduced, with 0 referring to the inertial and 1 to the corotating frame (compare Fig. 1). Inertial and corotating coordinates are related to each other via a passive rotation,

(4)$${}^0{\boldsymbol \omega} = {\bf R}\, {}^1{\boldsymbol \omega}$$

An efficient algorithm for the integration of (2) relies on the midpoint rule,

(5)$${\bf R}_{n + 1} = \exp\left(\Delta t\, {}^0{\bf \Omega}_{n + {1\over 2}} \right){\bf R}_n$$

with Δt = t _n+1 − t _n being the time step size. The update scheme (5) ensures that the orthogonality of the rotation matrix is conserved.Footnote ⁶ An obvious choice for the mid-point angular velocities is given by

(6)$${}^0{\boldsymbol \omega}_{n + {1\over 2}} = {1\over 2}\left({}^0{\boldsymbol \omega}_n + {}^0{\boldsymbol \omega}_{n + 1} \right) = {1\over 2}\left({\bf R}_n\, {}^1{\boldsymbol \omega}_n + {\bf R}_{n + 1}\, {}^1{\boldsymbol \omega}_{n + 1} \right)$$

where ¹ω_n and ¹ω_n+1 are the measured, corotational angular velocity components at the n-th and (n + 1)-th time step, respectively. Equation (6) delivers an implicit definition of components ⁰ω_n+1/2, which becomes apparent when inserting (5).

To derive additional advantage from the magnetometer data, the ansatz (6) is modified as follows. Use is made of the fact that the components of earth's magnetic field with respect to the global frame, ⁰B, are virtually constant in time,

(7)$${\bf R}\, {}^1{\bf B} = {}^0{\bf B} = {\rm const}$$

In particular, R_n+1 ¹B_n+1 = ⁰B. It is important to have in mind that ||¹B|| = ||⁰B|| = : B ₀.

Inserting (5), a non-linear system of equations constraining the three components ⁰ω_n+1/2 is obtained,

(8)$$\exp\left(\Delta t\, {}^0{\bf \Omega}_{n + {1\over 2}} \right){\bf R}_n\, {}^1{\bf B}_{n + 1} = {}^0{\bf B}$$

Note that the system (8) is under-determined. Together with the three equations (6), an over-determined system is obtained, from which the three components ${}^0{\boldsymbol \omega }^\ast _{n + {1}/{2}}$ can be determined via the least squares method,Footnote ⁷

(9)$${\bf F}\left({\boldsymbol \omega}^\ast \right)\equiv \left[\matrix{\; \, \Delta t\left\{\exp\left(\Delta t\, {\bf \Omega}^\ast \right){\bf R}_n\, {}^1{\boldsymbol \omega}_{n + 1} + {}^0{\boldsymbol \omega}_n - 2\, {\boldsymbol \omega}^\ast \right\}\cr \displaystyle {w_{n + 1}\over B_0}\left\{\exp\left(\Delta t\, {\bf \Omega}^\ast \right){\bf R}_n\, {}^1{\bf B}_{n + 1} - {}^0{\bf B} \right\}}\right] = {\bf 0}$$

Both systems, (6) and (8), have been scaled, such that they are free of any physical unit. In addition, a factor w _n has been introduced, which allows to weight the two subsystems relatively to each other. This turns out to be beneficial, since the reliability of the measured magnetic field is not constant throughout the measurement period. The respective reliability can be estimated from the deviation of the norm ||¹B_n|| from its expected value B ₀. On the basis of numerical experiments the following weighting function has been established,

(10)$$w_n = {\rm e}^{-\left(\,p\, {B_0-\Vert {}^1{\bf B}_n\Vert \over B_0} \right)^2}.$$

The parameter p is determined such that the noise of the obtained global magnetic field components is virtually the same for all three components ⁰B _x,y,z(t), see Figure 10, yielding a value of p = 5. For each time step n, the system (9) can be solved for the ω* by means of the Gauss–Newton method (Björck, Reference Björck1996). The rotation matrix is then updated via (5) with ${}^0{\boldsymbol \omega }_{n + {1\over 2}} = {\boldsymbol \omega }^\ast$.

With the help of the updated rotation matrix R_n+1, the measured quantities ¹B_n+1 (magnetic field), ¹ω_n+1 (angular velocity) and ¹a_n+1 (acceleration) are transformed to the global coordinate system via a passive rotation,

(11)$${}^0\{ {\bf B},\; {\boldsymbol \omega},\; {\bf a}\} _{n + 1} = {\bf R}_{n + 1}\, {}^1\{ {\bf B},\; {\boldsymbol \omega},\; {\bf a}\} _{n + 1}$$

To specify the initial orientation of the sensor unit, the global coordinate system $( {\bf e}^0_x,\; \, {\bf e}^0_y,\; \, {\bf e}^0_z)$ is defined according to the geographic East-North-Up (ENU) convention: The global z and y axes are upward vertical and geographic north, respectively. At rest, the IMU measures an acceleration vector ¹a₀ = −¹g, with g being earth's gravitational acceleration vector and a magnetic field ¹B₀, which is considered to originate solely from earth's magnetism. Correspondingly (Yun and others, Reference Yun, Bachmann and McGhee2008),

(12)$${\bf e}_x = {{}^1{\bf B}_0\times{}^1{\bf a}_0\over \Vert {\bf B}_0\times{\bf a}_0\Vert },\; \quad {\bf e}_y = {\bf e}_z\times{\bf e}_x,\; \quad {\bf e}_z = {{}^1{\bf a}_0\over \Vert {\bf a}_0\Vert }$$

The vector components ¹a₀ and ¹B₀ have been determined by averaging the measured values over a time period of 0.2 s immediately before the sensor unit has been entrained by the avalanche. To be more precise, the declination D of earth's magnetic field has to be taken into account, such that

(13)$${\bf e}^0_x = \cos D\, {\bf e}_x - \sin D\, {\bf e}_y,\; \quad {\bf e}^0_y = \sin D\, {\bf e}_x + \cos D\, {\bf e}_y,\; \quad {\bf e}^0_z = {\bf e}_z$$

A value of $D = 2.09705^\circ$ has been extracted according to Schnegg (Reference Schnegg1998). Correspondingly, the initial orientation is given by the rotation matrix

(14)$${\bf R}_0 = \left[\matrix{{\bf e}^0_x & {\bf e}^0_y & {\bf e}^0_z }\right]^T$$

which is used to determine ⁰B = R₀ ¹B₀ and to initialize the integration procedure.

4.1.2. Recovery of translational velocity and position

It is essential to distinguish between the ‘kinematic’ acceleration vector, here denoted as ${}^0\dot {\bf v}_n$, and the ‘dynamic’ acceleration vector ¹a_n, whose components are measured by the accelerometers. The former has to be understood as negative inertial force vector per mass, and the latter as guiding (or reaction) force per mass:Footnote ⁸ At rest, the kinematic acceleration is zero, ${}^0\dot {\bf v} = {\bf 0}$, whereas an ideal accelerometer triad is expected to measure the reaction to Earth's acceleration, ¹a = −¹g. In the course of an ideal ballistic motion, the kinematic acceleration equals Earth's acceleration, ${}^0\dot {\bf v} = {}^0{\bf g} = [ \, 0\; \; 0\; \; -g\, ] ^T$, whereas inertial forces and earth's attraction cancel each other, such that ¹a = 0. Generally,

(15)$${}^0\dot{\bf v}_n = {}^0{\bf g} + {\bf R}_n{}^1{\bf a}_n$$

with

(16)$${}^0{\bf g} = -{\bf R}_0\, {}^1{\bf a}_0 = \left[\matrix{0 & 0 & -\Vert {}^1{\bf a}_0\Vert }\right]^T$$

For a hypothetically perfect calibration of the accelerometers, one expects ||¹a₀|| to match the local value of Earth's gravitational acceleration (Schwartz and Lindau, Reference Schwartz and Lindau2003).

Equation (15) delivers the kinematic acceleration of the reference point, to which the measured accelerations are referred (see Fig. 1). The origin of the sensor coordinate system is arbitrarily chosen to be located at the geometric centre of the (almost) spherical housing. The respective coordinates of the reference point are

(17)$${}^1{\boldsymbol \varepsilon} = \left[\matrix{12 & -31 & -12 }\right]^T$$

with an estimated accuracy of ±1 mm (Fig. 1). Taking this eccentricity into account, the kinematic acceleration vector of the centre point is obtained by compensating for the rotational inertia term, $\dot {\boldsymbol \omega }\times {\boldsymbol \varepsilon }$, and a centripetal acceleration, ${\boldsymbol \omega }\times ( {\boldsymbol \omega }\times {\boldsymbol \varepsilon })$, see Wittenburg (Reference Wittenburg2008),

(18)$${}^0\dot{\bf v}_{n + 1}^{( C) } = {}^0{\bf g} + {\bf R}_{n + 1}\left[{}^1{\bf a}_{n + 1} - {}^1\dot{\boldsymbol \omega}_{n + 1}\times{}^1{\boldsymbol \varepsilon} - {}^1{\boldsymbol \omega}_{n + 1}\times( {}^1{\boldsymbol \omega}_{n + 1}\times{}^1{\boldsymbol \varepsilon}) \right]$$

The angular acceleration $\dot {\boldsymbol \omega }$ can be approximated by a symmetric finite difference,

(19)$${}^1\dot{\boldsymbol \omega}_{n + 1} \approx {{}^1{\boldsymbol \omega}_{n + 2} - {}^1{\boldsymbol \omega}_n\over t_{n + 2} - t_n}$$

An estimate of sensor unit velocities is obtained by integrating the kinematic accelerations (18) via the trapezoidal rule with ⁰v₀ = 0,

(20)$${}^0{\bf v}_{n + 1} = {}^0{\bf v}_n + \Delta t\, {{}^0\dot{\bf v}_n + {}^0\dot{\bf v}_{n + 1}\over 2}$$

and, analogously, for the positions,

(21)$${}^0{\bf x}_{n + 1} = {}^0{\bf x}_n + \Delta t\, {{}^0{\bf v}_n + {}^0{\bf v}_{n + 1}\over 2}$$

To understand the role of eccentricity and to quantitatively assess its effect on velocity and position results, it is helpful to imagine what the results of a hypothetically perfect tracking procedure (neither measurement nor integration errors) are expected to be: such procedure would deliver the trajectory of the centre point according to the translational motion superimposed by the spiralling path of the accelerometer reference point according to the rotational motion. Consequently, the trajectories of centre and reference point, respectively, would not deviate from each other by more than a distance of $\Vert {}^1{\boldsymbol \varepsilon }\Vert$. This gives rise to the assumption, that also under realistic circumstances, the influence of a small eccentricity on the trajectory is small, since the oscillating disturbances related to the extra terms in (18) will largely be suppressed by the smoothing effect of the numerical integration, as long as no parasitic aliasing occurs related to the specific sampling rate. This assumption must, of course, be verified on the basis of the experimental results, see the next section.

4.2.1. Slope coordinate system

It is a favourable circumstance that the relevant section of the terrain which covers the first 9 s of the sensor unit motion is virtually planar and can thus be approximated by an inclined plane with an inclination of $\alpha \approx 32^\circ$. The projection of the vertical axis ${\bf e}_z^0$ onto this plane will be called the slope secant. The name derives from the fact that it is determined from the TLS terrain model and passes through the starting point of the sensor unit and a corresponding point at the terrain surface in a distance of 120 m, which is approximately the travelling distance after 9 s. The (exactly straight) slope secant approximates the (almost straight) slop line, which is the projection of the vertical axis onto the actual terrain surface. Correspondingly, an additional Cartesian coordinate frame $( {\bf e}_x,\; \, {\bf e}^s_y,\; \, {\bf e}^s_z)$, which is particularly useful for the interpretation of the tracking results, can be defined in a straight-forward manner: the first basis vector ${\bf e}^s_x$ is aligned with the slope secant. The second basis vector ${\bf e}^s_y = {\bf e}_z^0\times {\bf e}^s_x$ is horizontal and tangent to the corresponding contour line. The third one, ${\bf e}^s_z = {\bf e}^s_x\times {\bf e}^s_y$ represents the direction normal to the terrain.Footnote ⁹ The so defined basis will be called the slope coordinate system (or frame) and its axis directions are referred to as longitudinal (${\bf e}^s_x$), lateral (${\bf e}^s_y$) and transverse (${\bf e}^s_z$). The rotation tensor transforming the inertial frame $( {\bf e}_x^0,\; \, {\bf e}_y^0,\; \, {\bf e}_z^0)$ into the slope frame is given by

(22)$${\bf Q}_s = \left[\matrix{{\bf e}^s_x & {\bf e}^s_y & {\bf e}^s_z }\right].$$

4.2.2. Rotational motion

The measurement range of the gyrometers, which has been chosen too small, has lead to a cut-off of a few angular velocity peaks somewhat above ±6 rad s${^{-1}}$ (which is, in fact, larger than the nominal range of $\pm 300^\circ$ s⁻¹). Luckily, whenever this has happened, not more than one out of the three angular velocity components has been affected by saturation. Therefore, the missing component can be recovered with the help of the magnetometer data. The corresponding procedure is described in detail in Appendix B. The results are shown in Figure 6 (left). The recovered curves provide a smooth continuation of the measured ones.

An estimate R_n for the angular orientation of the sensor unit at any time step t _n in the considered time interval is obtained by integrating the calibrated gyrometer data (Eqn (5)). For this purpose, the mid-point values of angular velocity components, ⁰ω_n+1/2, are calculated by solving the overdetermined, non-linear equation system (9) via the Gauß–Newton method. The integration procedure is initialized with the starting orientation given by (14). It shall be mentioned here (and will be discussed subsequently), that the recovered angular orientation relies not only on angular rates, but also on magnetometer data and on accelerometer data (via the initial orientation). The rotational motion so obtained can be visualized and interpreted favourably via the evolution of the rotation (or Euler) vector φ_n, which is implicitly defined by

(23)$${\bf R}_n{\bf R}_0^T = \exp{\bf \Phi}_n$$

Wittenburg (Reference Wittenburg2008) with Φ_n being the skew-symmetric matrix composed from the components of vector φ_n, compare Eqn (3). The matrix ${\bf R}_n{\bf R}_0^T$ refers to the orientation at time step t _n relatively to the initial orientation R₀. It shall be noticed that the components R_n and φ_n refer to the global, inertial coordinate system $( {\bf e}_x^0,\; \, {\bf e}_y^0,\; \, {\bf e}_z^0)$ as well as to the local, moving coordinate system $( {\bf e}_x^1,\; \, {\bf e}_y^1,\; \, {\bf e}_z^1)$ at the same time, i.e. ⁰R_n = ¹R_n and ⁰φ_n = ¹φ_n (Wittenburg, Reference Wittenburg2008).Footnote ¹⁰ For a straight-forward interpretation of the rotational motion it is, however, beneficial to express the rotation vector with respect to the slope frame $( {\bf e}_x^s,\; \, {\bf e}_y^s,\; \, {\bf e}_z^s)$, indicated through the left superscript symbol ‘s’,

(24)$${}^s{\boldsymbol \varphi}_n = {\bf Q}_s^T{\boldsymbol \varphi}_n$$

The result is drawn in Figure 7 and compared to the angular velocity components ${}^s{\boldsymbol \omega }_n = {\bf Q}_s^T{}^0{\boldsymbol \omega }_n$ with respect to the same frame. It can be concluded from these diagrams that in the course of the first 3.2 s the rotational motion is dominated by a fast rotation about the horizontal ${\bf e}^s_y$-axis at an average speed of $7.81\, {\rm rad}\, {\rm s} {^{-1}}\approx 477^\circ \, {\rm s}\, {^{-1}}\approx 1.24$ rounds per second. In this time period, the sensor unit performs a number of four full twits about the ${\bf e}^s_y$-axis. The first twist is the fastest one with an average speed of $12.36\, {\rm rad}\, s {^{-1}}\approx 708^\circ \, {\rm s} {^{-1}}\approx 1.97$ rounds per second and an instantaneous velocity peak of approximately $35\, {\rm rad}\, {\rm s} {^{-1}}\approx 2000^\circ \, {\rm s} {^{-1}}\approx 5.5$ rounds per second.Footnote ¹¹ For t > 3.2 s, the rotational motion continues at a lower speed and without a stable rotation axis. As mentioned before, the results for t > 9 s are corrupt due to a data gap and are, thus, not displayed here. A demonstrative visualization of the rotational motion is presented in Figure 8. For the visualization, the path r(t) drawn there is the trajectory of a point on a hypothetical sphere with an (arbitrarily chosen) radius of r = 5 m which is thought to move and rotate in the same way as the sensor unit does, i.e.

$${\bf r}( t_n) = {\bf x}_n + {\bf R}_n\left(r{\bf e}_z^0\right)$$

Figure 7.

Rotation vector (left) and angular velocity (right) components with respect to the slope frame.

Figure 8.

Red: Trajectory of the centre of the sensor housing, x(t), corresponding to the time interval [0,  9] s. Blue: Trajectory of a hypothetical tracing point in a constant distance of 5 m from the centre, r(t), for visualization of the rotational motion. Left: Axonometric view. Right: Top view, i.e. a normal projection of the trajectory onto a horizontal plane.

4.2.3. Magnetic field

It has been described in detail how the magnetometer data can be used to stabilize the integration of angular rates (Section 4.1.1) and to bridge short gaps with deficient gyrometer data (Appendix B). In addition, the measured magnetic flux density can give an evidence to assess the reliability of the recovered angular orientations. In this context, a proper calibration is crucial. At the absence of ferromagnetic disturbances it can be assumed that the vector of the magnetic flux density, and thus its components with respect to the inertial frame, ⁰B, are virtually constant in space and time. Consequently, the norm of the local components ||¹B|| is expected to be constant as well,

(25)$$\Vert {}^1{\bf B}( t) \Vert = \Vert {}^0{\bf B}( t) \Vert = {\rm const}$$

A sufficiently accurate fulfilment of Eqn (25) is thus a necessary condition for the measured values and, in particular, for the calibration to be reliable. A significant deviation from the expected value is used to define a quantitative reliability measure according to Eqn (10), which governs the effect of the magnetometer data on the angular orientations, compare Figure 9.

Figure 9.

Left: Vector norm of the local components of the magnetic flux density, ‖¹B‖, subjected to a calibration under laboratory conditions carried out some time after the avalanche experiment (cyan) and to an in situ calibration according to Appendix B (blue). The reliability measure (gray) is calculated from (10). Right: Local components of the magnetic flux density, ¹B_k, subjected to in situ calibration.

Figure 10.

Left: Inclination I(t _n) calculated from ⁰B(t _n) via Eqn (27b), its mean value $\bar I$ and the true value I ₀. Right: Global components of the magnetic flux density, ⁰B _x,y,z. The B components involved are subjected to in situ calibration, whereas the initial orientation relies on accelerometer data subjected to laboratory calibration.

A necessary condition for the reliability of the recovered orientations R(t) is given by the constancy of the global components of B,

(26)$${}^0{\bf B} = {\bf R}( t) {}^1{\bf B}( t) = {\rm const}$$

which is illustrated in Figure 10. It is also worth to discuss the direction of vector B specified via spherical (azimuthal and polar) angles, φ and $\vartheta$. Defining these angles with respect to the global ENU coordinate system $( {\bf e}^0_x,\; \, {\bf e}^0_y,\; \, {\bf e}^0_z)$, one has

(27)$$\varphi = \arctan{{}^0\!B_y\over {}^0\!B_x} = 90^\circ - D,\; \quad \vartheta = \arccos{{}^0\!B_z\over \Vert {\bf B}\Vert } = 90^\circ + I$$

with declination D and inclination I.

According to Schnegg (Reference Schnegg1998), at the time and location of the experiment, the corresponding values are $D_0 = 2.09705^\circ$ and $I_0 = 62.81433^\circ$. The declination has already been used for the specification of the starting orientation in (13) and can therefore not serve for any further verification. The inclination, on the other hand, has not been considered so far. Calculating the latter from the global magnetic field components via Eqn (27b), values are obtained which oscillate around a mean value of $\bar I \approx 61.19^\circ$ with a standard deviation of $1.3^\circ$ (see Fig. 10). Thus, a deviation of $I_0-\bar I\approx 1.62^\circ$ from the expected value I ₀ is observed. This result reveals a moderate error in the recovered orientations.

4.2.4. Translational motion

In Figure 11, at the left-hand side, the dynamic acceleration components with respect to the sensor frame, ¹a, are drawn. These are the values measured by the accelerometers and subjected to laboratory calibration. At the right-hand side, the kinematic acceleration components with respect to the inertial frame, ${}^0\dot {\bf v}$, Eqn (18), are plotted. It is emphasized that, in contrast to the ¹a, the ${}^0\dot {\bf v}$ rely on the whole set of measurement data, since the rotation matrix R is required for the coordinate transformation (11).

Figure 11.

Left: Dynamic acceleration components with respect to the local coordinate system. Right: Kinematic acceleration components with respect to the global coordinate system.

The clear representation of ballistic motion is considered as a certain evidence for the reliability of the measurements: a pronounced ballistic phase starts at t ≈ 0.85 s and is followed by several shorter ones, especially in the time period from t ≈ 7.6 s to t ≈ 9 s, during which a sort of saltational motion occurs (short ballistic phases interrupted by impacts; see Fig. 12). The ballistic phases are represented by the expected values of ${}^0\dot {\bf v} = [ 0\; 0\; -g] ^T$.

Figure 12.

Kinematic acceleration components related to the global frame: selected time intervals dominated by ballistic or saltational motion.

Knowing the evolution of the kinematic acceleration, the sensor unit trajectory can be recovered according to Section 4.1.2. Projections of the obtained trajectory and of the reference tube (Section 2) are depicted in Figure 13. The plots refer to the whole time span [0,  16] s in order to demonstrate the parasitic effect of the data gap at t ≈ 9 s. Correspondingly, the second third of the path (dashed lines) might be dismissed. Note that the reference domain applying to the top view (left diagram) is also drawn in Figure 2 and recall that the reference domain for the vertical section (right diagram) is aligned with the terrain surface (including the snow cover) and has an assumed vertical span of ±1 m.

Figure 13.

Sensor unit trajectory recovered from IMU data according to Section 4.1.2 and compared to the reference defined in Section 1. The heavy deviations from the reference in the second third of the path (dashed lines) are a result of the data gap at t ≈ 9 s and might be dismissed. Left: Top view, i.e. a normal projection onto a horizontal plane. The blue and red curves refer to an evaluation of (18) with and without consideration of eccentricity, respectively. Right: Vertical section, i.e. a normal projection onto a vertical plane passing through the slope line.

Figure 13 also demonstrates the influence of the eccentricity ${}^1{\boldsymbol \varepsilon }$ included in Eqn (18): the blue and red curves refer to an evaluation of (18) with and without consideration of eccentricity, i.e. with ${}^1{\boldsymbol \varepsilon }$ given by (17) and with ${}^1{\boldsymbol \varepsilon } = {\bf 0}$, respectively. The small extend of the deviation of the corresponding results indicates that the considerations of Section 4.1.2 concerning the effect of eccentricity are applicable. In particular, the sensitivity of the tracking results to a small change of the eccentricity vector is negligible. Accordingly, an eccentricity given by (17) is considered for all subsequent computations.

A comparison of the recovered trajectory with the reference reveals a somewhat unexpected behaviour: referring to the top view in Figure 13 (left), the numerical result matches the reference surprisingly well. A slight but progressive lateral movement is observed. The most likely explanation for this behaviour is a (quadratic) numerical drift, which is expected due to the twofold integration procedure leading to the position results. However, it cannot be excluded, that (to a small extent) the lateral movement also reflects a physical motion. In contrast, the vertical section in Figure 13 (right) reveals a significant deviation from the terrain surface, which is physically impossible. This deviation is more like an overall misalignment rather than an expected drift, since it is present even at the beginning of the movement and is virtually constant throughout the path (until the data gap occurs at t ≈ 9 s). This observation is supported by the dotted red line in the right diagram, which results from a rotation of the solid red line by an assumed angle of $6^\circ$ in the vertical section plane.

5.1.1. Accelerometer error model

The recalibration of the three-axis accelerometer based on in situ data relies on a common, linear error model, similar to the one which has already been used for the standard calibration procedure (Appendix A): it is assumed that the ‘true’ acceleration values ¹a_n are related to the measured ones, ${}^1\tilde {\bf a}_n$, via

(28)$${}^1{\bf a}_n = \left({\bf I} + {\bf C} \right){}^1\tilde{\bf a}_n + g{\bf b}_0 + g\tau_n{\bf b}_1$$

Here, the calibration matrix A = I + C is assumed to deviate only slightly from the identity matrix I. It allows for scaling errors, non-orthogonality and misalignment of the physical sensor axes relatively to the axes of the sensor coordinate system. The additive offset consists of a constant part (bias) g b₀ and a linear drift g b₁. Earth's acceleration g has been premultiplied to achieve calibration parameters which are independent of the choice of physical units. Due to the same reason, a dimensionless time parameter is used,

(29)$$\tau_n = { t_n - t_0 \over t_N - t_0}\in [ 0,\; 1] $$

with N being the index of the last considered time step. The recalibration also affects ⁰g (Eqn (16)), which becomes

(30)$${}^0{\bf g} = -{\bf R}_0\, {}^1{\bf a}_0 = \left[\matrix{0 & 0 & -\left\Vert \left({\bf I} + {\bf C} \right){}^1\tilde{\bf a}_0 + g{\bf b}_0 \right\Vert }\right]^T$$

For the proposed implementation, the ‘true’ values are adopted, i.e. ⁰g = [ 0 0  − g ]^T with g = 9.802 m/s⁻² for the local region.Footnote ¹² The kinematic acceleration vectors (Eqn (15)) now read

(31)$${}^0\dot{\bf v}_n = {}^0{\bf g} + {\bf R}_n\left[\left({\bf I} + {\bf C} \right){}^1\tilde{\bf a}_n - g{\bf b}_0 - g\tau_n{\bf b}_1 \right]$$

By rearranging the calibration parameters

(32)$${\bf C} = \left[\matrix{c_1 & c_4 & c_7\cr c_2 & c_5 & c_8\cr c_3 & c_6 & c_9 }\right],\; \quad {\bf b}_0 = \left[\matrix{b_1 \cr b_2 \cr b_3 }\right],\; \quad {\bf b}_1 = \left[\matrix{b_4 \cr b_5 \cr b_6 }\right]$$

the vectors

(33)$${\bf c} = [ \, c_1\; \dots\; c_9\, ] ^T,\; \quad {\bf b} = [ \, b_1\; \dots\; b_6\, ] ^T$$

are obtained, by which means Eqn (31) can be given a concise shape,

(34)$${}^0\dot{\bf v}_n = {}^0\dot{\tilde{\bf v}}_n + {\bf A}_n^{\prime\prime}{\bf b} + {\bf B}_n^{\prime\prime}{\bf c}$$

Therein,

(35)$${}^0\dot{\tilde{\bf v}}_n = {}^0{\bf g} + {\bf R}_n{}^1\tilde{\bf a}_n$$

are the kinematic acceleration vectors ahead of recalibration. The (3 × 6)-matrices

(36)$${\bf A}_n^{\prime\prime} = g\left[\matrix{{\bf R}_n & \tau_n{\bf R}_n }\right]$$

impart the correction according to bias and drift, whereas the (3 × 9)-matrices

(37)$${\bf B}_n^{\prime\prime} = \left[\matrix{{\bf b}^{( n) }_{11} & {\bf b}^{( n) }_{21} & {\bf b}^{( n) }_{31} & {\bf b}^{( n) }_{12} & {\bf b}^{( n) }_{22} & {\bf b}^{( n) }_{32} & {\bf b}^{( n) }_{13} & {\bf b}^{( n) }_{23} & {\bf b}^{( n) }_{33} }\right]$$

are related to the correction of scaling, non-orthogonality and misalignment. They are composed of columns

(38)$${\bf b}^{( n) }_{ij} = {\bf R}^{( n) }_{\, \colon \, , i}\, {}^1\tilde a^{( n) }_j$$

where ${\bf R}^{( n) }_{\, \colon \, , i}$ denotes the i-th column of the rotation matrix R_n and ${}^1\tilde a^{( n) }_j$ the j-th component of vector ${}^1\tilde {\bf a}_n$.

The improved kinematic accelerations (34) yield improved velocity and position vectors,

(39)$${}^0{\bf v}_n = \int_{t_0}^{t_n}{}^0\dot{\bf v}( t) \, {\rm d}t = {}^0\tilde{\bf v}_n + {\bf A}_n^\prime{\bf b} + {\bf B}_n^\prime{\bf c}$$

(40)$${}^0{\bf x}_n - {}^0{\bf x}_0 = \int_{t_0}^{t_n}{}^0{\bf v}( t) \, {\rm d}t = {}^0\tilde{\bf x}_n + {\bf A}_n{\bf b} + {\bf B}_n{\bf c} - {}^0{\bf x}_0$$

Again, ${}^0\tilde {\bf v}_n$ and ${}^0\tilde {\bf x}_n$ denote the results ahead of recalibration, and

(41)$${\bf A}_n^\prime = \int_{t_0}^{t_n}{\bf A}\!{}^{\prime\prime}( t) \, {\rm d}t,\; \quad {\bf B}_n^\prime = \int_{t_0}^{t_n}{\bf B}^{\prime\prime}( t) \, {\rm d}t$$

(42)$${\bf A}_n = \int_{t_0}^{t_n}{\bf A}\!{}^\prime( t) \, {\rm d}t,\; \quad {\bf B}_n = \int_{t_0}^{t_n}{\bf B}^\prime( t) \, {\rm d}t$$

The integration is again performed via the trapezoidal rule. Correspondingly, the integral operations in Eqns (39)–(42) have to be understood as

(43)$$\int_{t_0}^{t_n} f( t) \, {\rm d}t\, \colon = \sum_{m = 1}^n{1\over 2}\, ( f_{m-1} + f_m ) \, ( t_m-t_{m-1}) $$

5.1.2. Topography constraint

In fact, the sensor unit has moved close to the topographic surface (terrain). Consequently, the recovered trajectory is also supposed to develop close to this surface. This gives rise to postulate a soft constraint either at the velocity level,

(44)$${\sum_{n = 1}^N} \left({\bf k}^T{\bf v}_n\right)^2 \quad \rightarrow \quad {\rm Min!}$$

or at the position level,

(45)$$\sum_{n = 1}^N \left[{\bf k}^T( {\bf x}_n - {\bf x}_0 ) \right]^2\quad\rightarrow\quad{\rm Min!}$$

with k being the vector normal to the topographic surface. Equations (44) and (45) claim that the component normal to the terrain either of the velocity or of the relative position is minimum in a least squares sense. To formulate these constraint conditions, the surface normal vector k is required, which can, in principle, be obtained from the TLS data. In the present case, k is almost constant and can be approximated by ${\bf k} = {\bf e}^s_z$ (Section 4.2.1). This circumstance makes the implementation significantly easier. It is, however, not a mandatory requirement for the method to work properly. Inserting (39) and (40), the constraint conditions (44) and (45) become

(46)$$\left\Vert \tilde{\bf v}_\perp + {\bf A}^\prime_\perp{\bf b} + {\bf B}^\prime_\perp{\bf c} \right\Vert ^2\quad\rightarrow\quad {\rm Min!}$$

(47)$$\left\Vert \tilde{\bf x}_\perp + {\bf A}_\perp{\bf b} + {\bf B}_\perp{\bf c} \right\Vert ^2\quad\rightarrow\quad {\rm Min!}$$

Conditions (46) and (47) are equivalent to the overdetermined systems of equations,

(48)$$\left[\matrix{{\bf A}^\prime_\perp & {\bf B}^\prime_\perp }\right]\left[\matrix{{\bf b} \cr {\bf c} }\right] = -\tilde{\bf v}_\perp,\; \quad \left[\matrix{{\bf A}_\perp & {\bf B}_\perp }\right]\left[\matrix{{\bf b} \cr {\bf c} }\right] = -\tilde{\bf x}_\perp$$

respectively. For conciseness, the following (N × 1)-vectors,

(49)$$\tilde{\bf v}_\perp = \left[\, {\bf k}^T\tilde{\bf v}_1\ \dots\ {\bf k}^T\tilde{\bf v}_N\, \right]^T$$

(50)$$\tilde{\bf x}_\perp = \left[\, {\bf k}^T( \tilde{\bf x}_1 - {\bf x}_0 ) \ \dots\ {\bf k}^T( \tilde{\bf x}_1 - {\bf x}_0 ) \, \right]^T$$

as well as (N × 6)- and (N × 9)-matrices,

(51)$${\bf A}^{( \prime) }_\perp = \left[\matrix{{\bf k}^T{\bf A}^{( \prime) }_1 \cr \vdots \cr {\bf k}^T{\bf A}^{( \prime) }_N }\right],\; \quad {\bf B}^{( \prime) }_\perp = \left[\matrix{{\bf k}^T{\bf B}^{( \prime) }_1 \cr \vdots \cr {\bf k}^T{\bf B}^{( \prime) }_N }\right]$$

have been introduced. Each of Eqns (48a) and (48b) constitutes an overdetermined, linear system of N equations for the 6 + 9 = 15 unknowns contained in vectors b and c. As discussed at the beginning of Section 5, the rotational motion involved must display a certain degree of complexity, to ensure that matrices (51a) and (51b) have full rank (i.e. 6 and 9, respectively). Consequently, each of the two equation systems yields a unique solution, which will be denoted (b_x,  c_x) for the position level constraint and (b_v,  c_v) for the velocity level.

5.1.3. Regularization

Determining the calibration parameters is a sort of inverse problem. In fact, none of the solutions of systems (48) is feasible, since at least some of the calibration parameters obtained (their absolute values) are too large, i.e. in the order of 1. This contrasts the original intention to search for a small correction of pre-calibrated values. Accordingly, the minimization problem has to be reformulated to search for a compromise between minimizing the overall deviation from the topographic surface and minimizing the calibration parameters. To this end, the objective functions are enriched by regularization terms,

(52)$$f_v( {\bf b},\; {\bf c}) \, \colon = \left\Vert \tilde{\bf v}_\perp + {\bf A}^\prime_\perp{\bf b} + {\bf B}^\prime_\perp{\bf c} \right\Vert ^2 + \rho\, g T\left(\Vert {\bf b}\Vert ^2 + \Vert {\bf c}\Vert ^2 \right)\rightarrow {\rm Min!}$$

(53)$$f_x( {\bf b},\; {\bf c}) \, \colon = \left\Vert \tilde{\bf x}_\perp + {\bf A}_\perp{\bf b} + {\bf B}_\perp{\bf c} \right\Vert ^2 + \rho\, {g T^2\over 2}\left(\Vert {\bf b}\Vert ^2 + \Vert {\bf c}\Vert ^2 \right)\rightarrow {\rm Min!}$$

This approach is well known as Tychonov regularization (Strang, Reference Strang2007). The scalings involving Earth's acceleration g and the span T of the considered time interval have been introduced to provide the respective regularization term with the correct physical unit. The optimal choice for the regularization parameter is ρ = 1, which will be verified in the next section (compare Fig. 14): firstly, the smallest minimum of the objective function can be achieved for this choice. Secondly, the minimum is relatively insensitive to a variation of the parameter in a certain interval around 1 (see Fig. 14). Again, the minimization problem (52) is equivalent to an overdetermined system, namely to the system of N + 6 + 9 equations,

(54)$$\left[\matrix{{\bf A}^\prime_\perp & {\bf B}^\prime_\perp \cr \rho{\rm g} T\, {\bf I}_6 & {\bf O}_{6\times9} \cr {\bf O}_{9\times6} & \rho{\rm g} T\, {\bf I}_9 }\right]\left[\matrix{{\bf b} \cr {\bf c} }\right] = -\left[\matrix{\tilde{\bf v}_\perp \cr {\bf 0}_6 \cr {\bf 0}_9 }\right]$$

and analogously for (53),

(55)$$\left[\matrix{{\bf A}_\perp & {\bf B}_\perp \cr \rho\, {g T^2\over 2}\, {\bf I}_6 & {\bf O}_{6\times6} \cr {\bf O}_{9\times6} & \rho\, {g T^2\over 2}\, {\bf I}_9 }\right]\left[\matrix{{\bf b} \cr {\bf c} }\right] = -\left[\matrix{\tilde{\bf x}_\perp \cr {\bf 0}_6 \cr {\bf 0}_9 }\right]$$

Therein, I₆ is the (6 × 6)-identity matrix, O_9×6 the (9 × 6)-zero matrix, 0₆ = O_6×1, etc.

Figure 14.

Left: Local acceleration components before (red/dark green/blue) and after (magenta/green/cyan) recalibration. Right: Obtained minima of the objective functions according to the minimization problems (52) and (53), respectively, in dependence of the regularization parameter ρ. Note that the minimum objective function can be considered as a measure for the violation of the constraint condition.

5.2.1. Determination of calibration parameters

The parameter vectors b and c of the accelerometer error model result from the postulate that the constraint condition is fulfilled in a least squares sense either at position level or at velocity level (Eqns (44) and (45), respectively). To this end, the overdetermined systems of equation (54) or (55) have to be solved. The input quantities entering this procedure are: R_n, ${}^1\tilde {\bf a}_n$, ${}^0\tilde {\bf v}_n$ and ${}^0\tilde {\bf x}_n$. The tilde refers to values based on the laboratory calibration (prior to recalibration).

Applying the velocity constraint without considering the linear drift term involving b₁, one obtains

$${\bf b}_0^{( v) } = \left[\matrix{-0.086263 \cr 0.011683\cr - 0.006500 }\right],\; \quad {\bf C}^{( v) } = \left[\matrix{-0.008369 & 0.025738 & -0.04306\cr - 0.010927 & 0.005685 & 0.039751 \cr 0.029010 & -0.086067 & 0.061000 }\right]$$

whereas the position constraint yields

$${\bf b}_0^{( x) } = \left[\matrix{-0.061926 \cr 0.023639 \cr 0.001839 }\right],\; \quad {\bf C}^{( x) } = \left[\matrix{-0.004843 & -0.002775 & -0.027261\cr - 0.013881 & 0.028386 & 0.014610 \cr 0.015133 & -0.044731 & 0.056375 }\right]$$

From a physical point of view, the two constraint conditions are expected to yield similar results. In fact, the values of corresponding parameters differ by less that ±20%,

$$\delta{\bf b} = \left[\matrix{-0.158\cr - 0.078 \cr - 0.054 }\right],\; \quad \delta{\bf c} = \left[\matrix{-0.017 & 0.133 & -0.074 \cr 0.014 & -0.106 & 0.118 \cr 0.065 & -0.194 & 0.022 }\right]$$

with the relative difference being defined as $\delta {\bf b} = ( {\bf b}_0^{( v) }-{\bf b}_0^{( x) }) /( \Vert {\bf b}_0^{( v) }\Vert + \Vert {\bf b}_0^{( x) }\Vert )$ and analogously for δ c. On a first glance, the impact of the recalibration on the acceleration components is small (see Fig. 14, left). Its impact on the recovered trajectories, however, is significant, as will be explicated in detail in the next section. The recalibration also effects the determination of the initial orientation R₀ (Eqns (12)–(14)). This effect can be quantified by identifying the rotation vector δ φ^(x) related to the rotation matrix $\delta {\bf R} = {\bf R}_{0, recal}{\bf R}_0^T$,

$$\delta{\boldsymbol \varphi}^{( x) } = \left[\matrix{-2.724^\circ & 0.620^\circ & -1.002^\circ }\right]^T$$

Thus, the observed inclination of the magnetic field changes according to

$$I_{x} = \bar I - \varphi^{( x) }_x \approx 61.191^\circ + 2.724^\circ \approx 63.92^\circ$$

That is, the correction of the inclination tends in the right direction, though it is too large, such that $I_0 - I_x\approx -1.10^\circ$ (instead of $+ 1.62^\circ$ prior to recalibration). Now the absolute error goes below the standard deviation ($1.3^\circ$). The given numbers refer to the position constraint. For the velocity constraint, the angular deviation is even more pronounced,

$$\delta{\boldsymbol \varphi}^{( v) } = \left[\matrix{-3.516^\circ & 2.475^\circ & -4.602^\circ }\right]^T$$

such that the error of the observed inclination becomes $I_0 - I_v\approx -1.89^\circ$, which is even (absolutely) larger than the error prior to recalibration.

Concluding, the influence of the regularization parameter ρ entering the minimization problems (52) and (53), respectively, shall be investigated (see Fig. 14, right). For the velocity as well as for the position constraint, the smallest minimum of the objective function is achieved for ρ = 1. Moreover, the minimum is relatively insensitive to a variation of the parameter in a certain interval around 1. There is a pronounced interval of insensitivity for the position constraint ($5\cdot 10^{-2}\lesssim \rho \lesssim 5$). This is a necessary condition for the regularization approach to be valid. For the velocity constraint, this interval is significantly smaller, which makes the corresponding results less robust. Note that the sensitivities of the minimum objective function with respect to a variation of the regularization parameter are scalar measures for the sensitivities of the position and velocity functions, x(t) and v(t), respectively. Referring to the interval 5 ⋅ 10⁻² ≤ ρ ≤ 5, the variation of the minimum objective function is in the range of ±0.02. Correspondingly, the variation of the x(t) function is of the same order of magnitude ($\sim \pm 2\%$).

5.2.2. Constrained trajectory

The calibration parameters obtained from the two different types of constraint conditions (and an accelerometer error model (28) involving a constant bias but no drift) are subsequently used to recalibrate the accelerometer values, which, in turn, are applied to calculate new estimates of the sensor unit trajectory. In Figure 15, these results are compared to each other, to the original estimate of the trajectory according to standard calibration on the basis of laboratory data, and to the reference tube. According to the design of the method, the vertical projections of the constrained trajectories (right diagram in the figure) match the terrain line in an averaged sense. In contrast, dealing with the vertical projection (left diagram), the overall deviation of the constrained trajectories from the reference tube appears to become larger compared to their unconstrained counterpart. For a more detailed assessment of the new results it is beneficial to consider the position coordinates with respect to the slope coordinate system,

$${}^s x = {\bf e}^s_x\cdot\left({\bf x} - {\bf x}_0\right),\; \quad {}^s y = {\bf e}^s_y\cdot\left({\bf x} - {\bf x}_0\right),\; \quad {}^s z = {\bf e}^s_y\cdot\left({\bf x} - {\bf x}_0\right)$$

That is, ^sx (longitudinal coordinate) is the distance from the starting point of the sensor unit measured along the slope secant (Section 4.2.1); ^sy (lateral coordinate) and ^sz (transverse coordinate) are the normal distances from the slope secant measured in horizontal direction and in the slope normal direction, respectively. With other words, ^sx is approximately the travelling distance, ^sy and ^sz measure the deviation from the slope secant. Accordingly, Figure 16 shows the deviation of the estimated trajectories from the slope secant. The following conclusions can be drawn from the diagrams:

1. 1. The acceleration phase ranging from ^sx = 0 m to ^sx ≈ 1 m is accompanied by a small lateral movement in the order of 0.1 m.

2. 2. The ballistic motion phase at the beginning of the sensor unit movement is clearly represented in the trajectory plot, where it appears in the shape of a ballistic parabola ranging from ^sx ≈ 1 m to ^sx ≈ 21 m (right diagram in Fig. 16).

3. 3. The reliability of the unconstrained trajectory is poor, since the orientation of the ballistic parabola is not plausible. (A correct starting orientation of the trajectory is considered as a necessary condition for the reliability of a tracking procedure.)

4. 4. The starting orientations of the constraint trajectories are within an acceptable range.

5. 5. Each of the trajectories constrained in transverse direction suffers from a severe (basically quadratic) drift in the lateral direction (left diagram), which is to be expected, given the low sampling rate applied (compared to contemporary standards).

6. 6. The transverse deviation from the slope line (right diagram) is of the same order of magnitude as the estimated accuracy of the reference tube until it starts to diverge close to the end of the processed time interval (corresponding to a travelling distance of ≈120 m).Footnote ¹³

7. 7. Based on the trajectory results, a clear preference for the position or the velocity constraint cannot be deduced. The explications of the previous section, however, assign a higher reliability to the position constraint.

Figure 15.

Sensor unit trajectories recovered from IMU data subjected to laboratory calibration (red) and to in situ calibration based on velocity (cyan) and position constraint (blue), respectively. The projections of the reference tube are shown in grey. Left: Top view, i.e. a normal projection onto a horizontal plane. Right: Vertical section, i.e. a normal projection onto a vertical plane passing through the slope secant.

Figure 16.

Estimates of the sensor unit trajectory drawn in the inclined x-y plane (i.e. approximately the terrain surface) of the slope coordinate system (left) and in the vertical x-z plane (right). The vertical grey lines indicates the end of the ballistic flight path. The horizontal lines passing through the origin refer to the slope secant. Correspondingly, the left and right diagrams show the lateral and transverse deviation from the slope secant. Note that the reference tube is derived from the measured terrain model and is therefore not aligned exactly with the slope secant.

The calculations have been repeated, once without considering any bias, i.e. b₀ = b₀ = 0 and once with all terms included in the error model (28), i.e. considering constant bias and linear drift and are referred to as

1. 1. No bias: ${}^1{\bf a}_n = ( {\bf I} + {\bf C} ) {}^1\tilde {\bf a}_n$

2. 2. Constant bias: ${}^1{\bf a}_n = ( {\bf I} + {\bf C} ) {}^1\tilde {\bf a}_n + g{\bf b}_0$

3. 3. Linear drift: ${}^1{\bf a}_n = ( {\bf I} + {\bf C} ) {}^1\tilde {\bf a}_n + g{\bf b}_0 + g\tau _n{\bf b}_1$

The corresponding tracking results are summarized in Figure 17. It turns out that models (2) and (3) yield virtually identical results, whereas model (1) does not work properly. This observation holds for the velocity as well as for the position constraint. Clearly, approach (2) is to be preferred since approach (3) relies on a larger number of parameters without delivering better results. Thus, all position and velocity results presented in this section without an explicit reference to the error model rely on approach (2).

Figure 17.

Comparison of error models (1)–(3) for the velocity (top) and position constraint (bottom).

5.2.3. Translational velocities

The recovered translational velocity components are plotted in Figure 18. In the left diagram, the global components obtained before (cyan, green, magenta refer to global z, x, y components, respectively) and after recalibration (red, darkgreen, blue) are compared to each other. Apparently, the vertical velocity component ⁰v _z is effected most by the recalibration. This is what has to be expected according to the nature of the constraint condition. The divergence of the results according to position and velocity constraints is less significant (solid vs. dotted lines). There is also a certain deviation in the evolution of the ⁰v _y components, which reflects the drift observed in the horizontal projection of the trajectory and is thus considered as a result of accumulated noise with no physical significance. For a further interpretation of the velocity results it is again useful to consider the components with respect to the slope coordinate system (Section 4.2.1), i.e. ${}^s{\bf v} = {\bf Q}_s^T{}^0{\bf v}$, which are plotted in the right diagram of Figure 18. Correspondingly, ^sv _x refers to the direction of the slope secant (longitudinal), ^sv _y, to the horizontal (lateral) direction and ⁰v _z to the direction normal to the terrain (transverse). As expected, the ^sv _y values are relatively small. However, instead of simply oscillating around zero, they reveal the linear drift mentioned before. The motion phases identified within the qualitative IMU data analysis are reflected by the velocity evolution:

1. 1. Acceleration phase from t ≈ 0.25 s to t ≈ 0.80 s: a strong and progressive increase of the ^sv _x component occurs (Fig. 18, right).

2. 2. Ballistic motion from t ≈ 0.80 s to t ≈ 1.67 s: horizontal components ⁰v _x and ⁰v _y are basically constant, while the vertical component −⁰v _z is linearly increasing according to Earth's acceleration (Fig. 18, left). The evolution of the transverse component ^sv _z indicates that a pronounced motion normal to the terrain surface occurs that changes from outward (^sv _z > 0) to inward (^sv _z < 0) orientation (Fig. 18, right).

3. 3. Impact at t = 1.67 s: a sudden deceleration occurs in longitudinal (${\bf e}^s_x$) as well as in transverse direction (${\bf e}^s_z$), see the right diagram of Figure 18.

4. 4. ‘Steady flow’ motion starting at t ≥ 1.67 s: the motion is dominated by the longitudinal component ^sv _x. It is slowly decelerating and is superimposed by a basically random jerking (Fig. 18, right).

Figure 18.

Translational velocity components vs. time. Global components (left) and components with respect to the slope frame (right). Comparison of results before and after recalibration. Solid lines refer to the position constraint, dotted lines to the velocity constraint.

On the basis of supplementary Doppler radar measurements, the motion of the avalanche front has been determined. The observed signal is shown in the diagram in Figure 3, from which the range-time curve of the avalanche front can be deduced (as a step function), see also the right diagram of Figure 19. Accordingly, one single velocity value for each of the subsequent range gates is determined. The span of the recovered sensor unit trajectory covers five entire range gates, which are represented by black framed error boxes in the left diagram of Figure 19. The length of each box corresponds to the length of the range gate (≈25 m, each) and the width is given by the estimated accuracy of the velocity evaluation (≈±1.5 m s⁻¹). For the representation in the diagram, position and length of each range gate as well as the measured velocity have been projected onto the slope secant.

Figure 19.

Translational velocity components with respect to the slope frame vs. travelling distance ^s x (left) and the corresponding path-time diagram (right). Both diagrams provide a comparison with the Doppler radar results (gray domains): velocities (left) and range-time diagram of the avalanche front (right) projected on the terrain. The horizontal span of the error boxes and bars refers to the length of each range gate (≈25 m). The vertical span of the error boxes refers to the uncertainty of the velocity evaluation (≈±1.5 m s⁻¹). Note that close to a travelling distance of 125 m (according to t ≈ 9 s) the data corruption occurs.

From the velocity evolution displayed in the left diagram of Figure 19, it can be seen that the avalanche is in a decelerating phase already before the sensor unit is entrained. It takes the sensor unit 10–20 m to reach a velocity, which is similar to the one of the avalanche front (≈20 m s${^{-1}}$), meaning that the sensor unit is necessarily several metres (definitely less than a 25 m range gate) behind the front at that instant. Subsequently, the sensor unit appears to move roughly at the same decelerating speed as the avalanche front. After a travelling distance of about 50 m, the sensor unit velocity seems to decrease less than the frontal velocity, i.e. the sensor unit moves towards the front, which it appears to reach after a travelling distance of 70–100 m. From then on, the sensor unit is expected to remain close to the avalanche front for a while. This plausible scenario fits well with the data presented in both diagrams of Figure 19. The comparison of sensor unit and frontal velocities on the same range axis is justified by the observation that the maximum distance between the sensor unit and the avalanche front is significantly smaller than the length of a range gate. In addition, the diagram supports the conclusion that the tracking procedure prior to recalibration overestimates the longitudinal velocity significantly, whereas the (‘constrained’) velocity according to in situ calibration better fits to the radar data. After a travelling distance of about 120 m (just before the onset of data corruption), the constrained velocity drops below the radar velocity. This behaviour is assessed to be unphysical and to originate from the escalating lateral drift occurring close to the end of the processed time interval.