Track frame approach for heading and attitude estimation in operating railways using on-board MEMS sensor and encoder

Measurement 184 (2021) 109898 Available online 29 July 2021 0263-2241/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Measurement journal homepage: www.elsevier.com/locate/measurement Track frame approach for heading and attitude estimation in operating railways using on-board MEMS sensor and encoder Eduardo Brialesa,∗, Pedro Urdab, José L. Escalonab aVirtualmechanics, S.L., Spain bDepartment of Mechanical and Manufacturing Engineering, University of Seville, Spain ARTICLE INFO Keywords: 3D orientation estimation Fusion algorithm Inertial sensors Experimental setup Railway mechanics Industrial environment ABSTRACT In this work, the orientation of a railway is estimated with a novel methodology based on multibody system kinematics using the railway-specific track frame. The proposed method improves the prediction model by considering the translational accelerations due to the track negotiation. To this end, the forward velocity of the vehicle, measured with an encoder, and the design geometry of the track are used. This algorithm has been tested on an operational underground light-metro railway with quite good results compared with other data fusion algorithms embedded in commercial Inertial Measurements Units (IMU) that contains no information about the real application whatsoever. 1. Introduction The estimation of the attitude and heading1of a rigid body in the 3D space from sensors signals has been a concern of main interest for several purposes in the recent years. While high-precision inertial sensors have been historically used for orientation estimation in vehicles navigation, the development of low-cost micro electromechanical systems (MEMS) technology has allowed the use of MEMS devices for navigation applications. These sensors that estimate orientation through acceleration and angular velocity are usually referred to as inertial measurement units (IMU) or 9 degrees-of-freedom sensors (9DOF or NDOF) if they also incorporate a magnetometer. Some others sensors used for these applications are Global Navigation Satellite Systems (GNSS), Computer Vision (CV) and encoders to determine forward velocity of the moving wheeled objects. In literature, examples of different estimation algorithms can be found for bio-logging [1], agriculture machines [2–4] and navigation of aerial [5,6], maritime [7] and ground vehicles [8–13] or aerial, maritime and ground vehicles altogether [14]. While these articles differ in the methodology used, they are all based on the integration of the output signals, usually the angular rate, and its correction with gravity and/or magnetic field (MARG) signals. The correction is necessary because the sole integration of the data leads to inaccuracies that tend to be accumulative in time. Different techniques can be applied to correct these inaccuracies. The source of the errors in low-cost IMU include non-orthogonalities in ∗Corresponding author. E-mail address: [email protected] (E. Briales). 1Attitude and heading reference system (AHRS) describes the orientation (or posing) of an object in space. This terminology is preferred as it groups gravity-affected angles, the roll and pitch, in attitude, while leaving the gravity-invariant angle, the yaw, for the heading. the local axis, non-linearities and scale factors in the input–output transfer function within the sensor, random biases and noise [9,15,16]. Moreover, 9DOF MEMS also incorporate the magnetometer sensor that is subjected to magnetic noise (soft and hard iron distortions) that can be high for some industrial or underground environments turning this sensor unsuitable to estimate posing. These effects are inherent to the construction or physical properties of the MEMS and cannot be completely removed even after calibration [17], so the best options to increase precision are: to compute better algorithms that are robust to these errors or improve the physical model used for the prediction of the orientation estimation. Among the algorithms used to reduce orientation error estimated from an IMU rigidly attached to a body, one can find those based on estimators, state observers, numerical integration, filtering or optimization. These methodologies can all be augmented with reliability tests for the input measurements before using the signals. They all fuse the information of the magnitudes obtained by the sensors but with different objectives, such as minimizing a cost function (optimization), approximating the behaviour of a certain model (state observer) or exploiting different sensor strengths (complementary filter). These different methods aim to enhance the capabilities of the algorithms rather than improving the physical model of the object orientation. A brief description of the capabilities of these different algorithms is given below. https://doi.org/10.1016/j.measurement.2021.109898 Received 31 January 2021; Received in revised form 12 July 2021; Accepted 14 July 2021 Measurement 184 (2021) 109898 2 E. Briales et al. To enhance the capabilities of the algorithms there are different technical tools. Estimators aim to define a time-invariant quantity based on the measurements that tends to the real value of the orientation and is stable to noise with the passing time. State observers [7,18–20] aim to determine the state of the system, that may include the sensors biases. Some authors also incorporate some parameters of the models as states of the augmented system, such as the velocity of the object or the IMU errors parameters [9]. The fusion algorithms are based on combining the information of two or more sensors measuring different variables of the same system. As the noise is usually treated as a random variable, its effect in the resulting signal can be cancel out combining different sensors, what makes fusion algorithms also filters. Some of the more common filters used in literature for estimating orientation are Linear Kalman Filter [8], Non-linear Kalman Filter (such as Extended [10,12], Unscented [21] , Cubature or High Degree Kalman Filter [22]), Complementary Filter [3,23,24] or Unknown Input Filter [25]. These algorithms usually reach high-precision results that are used as reference values for checking new algorithms proposed in the literature. However, they tend to be computational costly because linear systems must be solved and their accuracy it not sufficient when the system model conditions are not met. Therefore, new more-efficient or ad-hoc algorithms are continually proposed. Finally, the optimization algorithms seek minimize – in one or several steps – the error function between the real orientation (unknown) and the estimations based on measured data. This optimization can be approached with any kind of minimization algorithm such as the Levenberg–Marquardt Algorithm (LMA) [13], Gauss–Newton Algorithm (GNA), Gradient Descent Algorithm (GDA) [26–28], Control Loop [6,29] and lately, some machine learning algorithms such as Least-Squares Support Vector Machine (LSSVM) and non-linear autoregressive with exogenous input (NARX) [13], Kalman filter based on Recurrent Neural Network (RNN) [30] or Fuzzy logic [31]. Criteria used to select among these available options are accuracy, reliability, computational efficiency, noise and drift stability, and robustness in the region of local minimum. To enhance these algorithm accuracy, there are some articles that defines reliability tests for the signals to be used before introducing them into the equations [11,32]. For example, the accelerometer only considered reliable when the acceleration of the movement is small compared to gravity. In addition to validation with experiments or simulations, some authors also include mathematical proofs of the stability of the proposed algorithms to noise using Lyapunov functions [1,6,18,33]. Additionally, to improve the physical model rather than the algorithms it is mandatory to focus on the application at hand, that in this work is the measurement of the attitude and heading of ground vehicles. Estimation algorithms have been tested for several ground vehicles obtaining measurements of the heading and attitude with promising results in different works. These works mainly focus on improving the prediction model including more information about the kinematics and dynamics of the motion or adding more sensors to the system. In many cases the use of Global Navigation Satellite Systems (GNSS) helps improving the orientation estimation such as in the works of [4,8,13,14,34,35] to cite some. However, GNSS cannot be used in applications such as light-metro underground vehicles with poor satellite connectivity. Some other authors make use of computer vision systems or laser sensors to add some information about attitude angles [36–38] that relies on periodic corrections based on features of interest. In the classic work for a ground robot, Vaganay and Aldon [10] use a combination of gravity and translational acceleration for the dynamic estimation of the rotation quaternion. In this case, the translational acceleration is estimated from rigid body kinematics equations and the result is extracted from the accelerometer output. Then, the accelerometer and gyroscope measurements are fused via a Kalman filter. Klier et al. [39] published another classic work that is frequently referenced for attitude estimation of ground vehicles. The authors formulate the kinematic problem as an Ordinary Differential Equation (ODE) using the translational velocity of the inertial sensor. They consider Euler angles and their derivatives as variables that are then implemented in a Kalman Filter. That work also implements some practical constrains for the velocity considering the usual ground vehicle motion such as small lateral velocity while driving straight or very small vertical velocity. This formulation of the problem has been widely replicated, for example, in the work of Scholte et al. [9], in which the ODE system implemented as an observer for state (sensor velocity and Euler angles) and parameter estimation (sensor biases and other geometric and dynamic parameters). In this case, the prediction model is augmented with force and moment balance equations of the vehicle, that are included as constrains. Marco et al. [16] use Klier’s model as baseline and the system is augmented with new geometric parameters, such as the slide-slip, and the dynamic balance equations once again. Similar approaches can be found in [34,40,41]. These authors also discuss about the observability of the model under certain navigation conditions (driving straight). Yuan et al.’s work [42] has a particularly interesting approach considering the effect of translational acceleration through the uncertainties of the observation model. However, it has been found that, although translational acceleration is one of the main reasons for the inaccuracy of attitude estimation methods, it has not been studied deep enough [43]. In the particular case of railways, that ride along a guideway with known geometry, this issue has not been properly addressed. In this work, a better prediction model for railways attitude and heading estimation is presented. It can also be applied to any other ground vehicles in which there is information about the guideway that follow. Narrowing the scope of application of the algorithm allows the implementation of a specific methodology based on railway kinematics which grants very accurate results in the estimation. Ultimately, this method is applied to obtain the orientation of the wheelset of a railway during a ride over an irregular track using solely an IMU and an encoder. In order to accomplish this task, the kinematic model of the vehicle moving along a track with known design geometry must be included in the approach. This kinematic model is implemented with the use of an intermediate frame (the so-called track frame) to isolate the gravity component in the accelerometer output measurements. Two different algorithms are developed. First, an implementation of Madgwick’s algorithm [26,27] based on a GDA that aims to obtain high accuracy while keeping a good efficiency. Second, an Extended Kalman Filter (EKF) that uses the gyro signals as inputs and the accelerometer signals as measurements [44,45]. This paper is organized as follows. Section 2shows the motivation for this work after finding inaccurate results in the orientation obtained from commercial IMU sensors and other Kalman Filter approaches based on inertial measurements with general purpose models. Section 3 introduces the different frames involved in the multibody system kinematics approach (global, track and body frames) as well as the Euler angles used. Proposed estimation model is presented in Section 4and then implemented in the GDA and EKF in Section 5. Section 6presents the equipment used for the data acquisition, the characteristics of the different experiments made and the data collected in them that are used as inputs of the algorithms. In Section 7the different results obtained with this new method are presented, implemented using different methods and compared. In Section 8some considerations about the experiments and hypotheses made are discussed and finally in Section 9 the main achievements of this method are summarized and the main conclusions are highlighted. 2. Motivation The interest on 3D orientation estimation emerges from a wider research in a railways application that aims to measure the irregular geometry of the track using solely inertial sensors, a computer vision system and an encoder, based on a kinematic multibody approach [46, 47]. One of the most influencing variables in the calculation of the irregular geometry is the rotation of the vehicle bodies with respect to Measurement 184 (2021) 109898 3 E. Briales et al. Fig. 1. Euler angles values obtained with commercial IMU LORD MicroStrain AHRS 3DM-GX4®-25. the track. Furthermore, the roll angle, 𝜑, is of utmost importance and it must be measured with the highest accuracy. By comparison, pitch angle has less influence in the final results and heading angle is not necessary for the geometry measurement. Unfortunately, in the area of experimental multibody systems, some of the most difficult magnitudes to be measured are the 3D orientation parameters. The inherent non-linear nature of rotation plus the high sensitivity of the Euler angles due to small changes in the measured acceleration make them difficult to be estimated. It is therefore hard to obtain values with enough precision given the aforementioned MEMS errors that lead to the integration-drift. The experience of the authors, after some field experiments on operational vehicles riding over irregular tracks, is that the values of the Euler angles provided by a high-quality commercial IMU2are not even close to the required accuracy for track geometry measurement. Fig. 1 shows the initial results granted by the mentioned commercial IMU and the track design angles. The solid red line shows the obtained Euler angles values while the dashed black line represent track design angles. Due to the IMU to track relative motion and due to the track irregularity, the IMU angles and the track design angles do not have to coincide, but the difference is not more than 1◦or 2◦. As it can be seen, the results for the three Euler angles used to characterize the rotation are far from the track design values in some sections. The angles accuracy obtained with the commercial IMU is worse than the specified in the sensor datasheet (±2◦for dynamic applications) reaching at some points 2.5◦of error for roll, 5◦in transient for pitch and more than 100◦of error for yaw. The IMU was programmed to use the magnetometer in order to include a reference for the heading. However, the presence of metal structures in the surroundings of the sensor may have caused high magnetic distortion [48–50] and it could have affected the poor output of the algorithm, mainly to the yaw angle. Some of the issues found in these values provided by the IMU are listed below: •For the roll (𝜑), the more noteworthy issue is the failure to properly measure the cant angles at curves. It is also striking the overshooting effect after the different transitions. •For the pitch (𝜃), it can be seen that there is a transient at the beginning of the signal probably due to the forward acceleration of the vehicle. Afterwards there are some issues in the signal at the points in which the other angles vary abruptly, for example at 16.35 km. 2LORD MicroStrain AHRS 3DM-GX4®-25. •For the yaw (𝜓), the reconstruction of the signal is completely inaccurate, failing to reconstruct even the trend of the reference magnitude. This could have been caused by different factors such as inadequate Euler angles sequence, linearization issues within the sensor, not-considered high translational accelerations, or magnetic distortions. When searching for this odd behaviour of the algorithm for the heading angle, little research was found for specific purpose applications in the railway sector or even for ground vehicles. The most likely cause was considered the magnetic distortions that may appear in industrial environments. Furthermore, one of the main hypothesis usually made for 3D orientation estimation is that the signals captured by the accelerometer are due mainly to the gravity. As the value of gravity is certainly high for a great deal of applications compared with the sensor accelerations, this hypothesis is usually valid. However, it is not the case for railways in which high frequency vibrations from the track irregularities and important lateral accelerations experimented by the vehicle when moving on curves may cause that the hypothesis conditions are not met. Before proposing alternative methods, it is helpful to replicate the usual implementation of the fusion algorithm supposedly based on Extended Kalman Filter (EKF) that might be used in different commercial IMU. Therefore, the orientation was recalculated implementing an EKF based on the works of Sabatini [44,45] in which the system model is based on the relation between the derivative of the angles and the angular velocities. The accelerometer signals are included afterwards in the update step, considering again that the acceleration signals are mainly determined by gravity. The obtained results are shown in Fig. 2. While the yaw calculation has greatly improved and the signal is practically the same than the reference (except for the small drift), the results for roll and pitch suffer from the same issues mentioned before. The issues for roll angle in the curves occur because translational acceleration was not considered. In this regard, there are some authors that include a certain information about this acceleration for general purpose applications. A particularly interesting work is Yuan et al. [42], in which this information is included as the uncertainty of the observation model in a EKF implementation. The results obtained with this approach are shown in Fig. 3. It can be seen that the behaviour during curves greatly improves for the roll angle, but the overall result is worse for pitch and heading angles. Moreover, a particular ripple appears in all the signals that is most likely due to the increased uncertainty for the measurement: even when the signal es filtered, there are mechanical vibrations that Measurement 184 (2021) 109898 4 E. Briales et al. Fig. 2. Euler angles obtained from commercial IMU LORD MicroStrain AHRS 3DM-GX4®-25 and comparison with implemented EKF. Fig. 3. Euler angles obtained with Yuan et al.’s model implementation in an EKF. are mixed with the acceleration values and increasing the uncertainty ripples the final result. This effect makes the Euler angles inaccurate mainly in straight sections of the track. These inadequate results force the pursuit of better approaches to estimate the 3D orientation of a railway travelling on a track. And this approach starts with a deeper understanding of the railway kinematics. 3. Kinematics In the multibody kinematics model used in this paper, the coordinates used to describe the motion of the bodies provide the position and orientation of the body frame (BF) with respect to an inertial global frame (GF). The model used in this work is based on that presented by the authors in [51], also used in [52]. As stated in that paper, for ground vehicles that follows a guideway (a track for trains) an intermediate frame that is neither inertial nor body-fixed can be defined for a convenient description of the motion of the bodies. This frame advances along the irregularity-free track with the same forward velocity as the vehicle. This intermediate frame is usually called track frame (TF) in railway dynamics (see Fig. 4). The TF is a function field of variable srather than a frame associated with a solid. Therefore, a particular TF can be defined for each BF, called body track frame, considering the TF whose 𝑋coordinate coincides with BF’s. The use of the GF is still necessary to obtain the equations of motion (EOM). The position and orientation of the TF with respect to the GF is a function of the track geometry and is usually obtained as a function of an arc-length coordinate. In the track frame formulation, body coordinates describe the position and orientation of the BF with respect to the TF. These coordinates are called here track coordinates. The main advantages of using the TF are: 1. Bounded magnitudes: The magnitude of the coordinates of the bodies remains bounded in the TF while they may grow with no limit when referred to the GF. 2. Accuracy: As a consequence of the previous point, the accuracy in the description of the relative motion of the vehicle with respect to the track is kept, even when the value of the absolute coordinates is large. Measurement 184 (2021) 109898 5 E. Briales et al. Fig. 4. Different frames involved in this kinematic approach: Global Frame,Track Frame,Body Frame. 3. Linearization: The formulation of kinematics equations in the TF allows the use of the small angles assumption that facilitates the linearization of trigonometric functions that describe the orientation. This allows to simplify the expressions and, therefore, the calculations. 4. Easier formulation of vehicle motion: Constrained motion equations of the track frame and vehicle bodies with respect to the track, like the forward velocity constrain, are easily imposed using track coordinates. On the other hand, the use of track coordinates shows the following drawbacks: 1. More complicated overall formulation: The kinematic description of the bodies is much more complicated because it involves the track design geometry. The calculation of velocities and acceleration of the bodies require the calculation of partial derivatives of the track geometry with respect to the TF arc-length coordinate. These expressions enter in the vehicle equations of motion becoming computationally more expensive. Moreover, the results usually rely on the accuracy of the track geometry data involved that sometimes cannot be guaranteed. 2. Application-bounded formulation: The inclusion of a railwayrelated frame, as the TF, hinders the use of general purpose multibody third-parties software. 3.1. Vectors and frames This section explains in detail the kinematic description used in the multibody track frame formulation used in this work for the modelling of the railway vehicles motion. In this work, the following symbols are used to express position vectors in different frames: 1.  𝑅is a position vector with respect to the GF. 2. 𝑟 is a position vector with respect to the TF. 3. 𝑢 is a position vector with respect to the BF. On the other hand, the following symbols are used to express the vector components given in the different frames: 1. Bold symbols without diacritic, like 𝐯, mean the 3 × 1 column matrix that contains the components of vector 𝑣 in the GF. 2. Bold symbols with ‘bar’ diacritic, like  𝐯, mean the 3×1column matrix that contains the components of vector 𝑣 in the TF 3. Bold symbols with ‘hat’ diacritic, like  𝐯, mean the 3 × 1 column matrix that contains the components of vector 𝑣 in the BF 3.2. Kinematics of the TF The TF [𝑂𝑡(𝑋𝑡𝑌𝑡𝑍𝑡)] is assumed to follow the track keeping axis 𝑋𝑡tangent to the track centreline. Axis 𝑌𝑡connects the two railheads centres. Therefore, the TF used in this work is not a Frenet frame of the track centreline since the chosen 𝑌axis is not parallel to the osculating plane. The track centreline is defined assuming an ideal reference geometry with absence of irregularities. The arc-length travelled by the TF along the projection in the horizontal XY plane of the ideal track centreline is denoted as 𝑠𝑡. This 𝑠𝑡coordinate does not follow the spacial 3D curve but its horizontal projection on the global plane ⟨𝑋𝑌 ⟩. The coordinates that describe the position of the origin and orientation of the TF with respect to the GF are: 𝐪𝑡=[𝑥𝑡𝑦𝑡𝑧𝑡𝜑𝑡𝜃𝑡𝜓𝑡]𝑇(1) where 𝐑𝑡= [𝑥𝑡𝑦𝑡𝑧𝑡]𝑇is the position vector of the origin of the TF, 𝑂𝑡, with respect to the GF and 𝚽𝑡= [𝜑𝑡𝜃𝑡𝜓𝑡]is a set of three Euler angles that describe the orientation of the TF with respect to the GF. The coordinates of the TF are known functions of the track centreline geometry, this is, 𝐪𝑡=𝐪𝑡(𝑠𝑡). In computational railway dynamics, these functions are implemented in the so-called track-preprocessor that uses as input the track design geometry parameters to build the functions 𝐪𝑡(𝑠). Therefore, the set of coordinates shown in Eq. (1) cannot be considered part of the system generalized coordinates. The arc-length coordinate, 𝑠𝑡, is in fact the only generalized coordinate that shows the gross forward motion of the vehicle along the track. The Euler angles sequence used in this work is the YPR, that stands for yaw,pitch,roll, the rotations about the 𝑍,𝑌′and 𝑋′′ axes, respectively. This sequence of angles results in the following rotation matrix: 𝐂=⎡ ⎢ ⎢ ⎢ ⎣ cos 𝜓cos 𝜃− sin 𝜓cos 𝜑+ sin 𝜃sin 𝜑cos 𝜓sin 𝜓sin 𝜑+ sin 𝜃cos 𝜓cos 𝜑 sin 𝜓cos 𝜃sin 𝜓sin 𝜃sin 𝜑+ cos 𝜓cos 𝜑sin 𝜓sin 𝜃cos 𝜑− sin 𝜑cos 𝜓 − sin 𝜃sin 𝜑cos 𝜃cos 𝜃cos 𝜑 ⎤ ⎥ ⎥ ⎥ ⎦ (2) Measurement 184 (2021) 109898 6 E. Briales et al. This Euler angles sequence is usually implemented for ground vehicles applications, as it leaves the gimbal lock singularity in an impossible rotation. The use of Euler angles instead of quaternions has some benefits such as better understanding the physics of the motion involved in the estimation of the orientation, as well as simplifying the linearization of the kinematics equations that includes the TF. 3.3. Velocities and accelerations of the TF The absolute velocity and acceleration of the TF is given by:   𝐑𝑡=⎡ ⎢ ⎢ ⎣ 𝑉 0 0 ⎤ ⎥ ⎥ ⎦ (3)   𝐑𝑡=⎡ ⎢ ⎢ ⎣  𝑉 𝜌ℎ𝑉2 −𝜌𝑣𝑉2⎤ ⎥ ⎥ ⎦ (4) And the absolute angular velocity and acceleration of the TF are given by:  𝛚𝑡=⎡ ⎢ ⎢ ⎣ 𝜌𝑡𝑤𝑉 𝜌𝑣𝑉 𝜌ℎ𝑉⎤ ⎥ ⎥ ⎦ (5)  𝛂𝑡=⎡ ⎢ ⎢ ⎣ 𝜌𝑡𝑤  𝑉 𝜌𝑣 𝑉 𝜌ℎ 𝑉+𝜌′ ℎ𝑉2⎤ ⎥ ⎥ ⎦ (6) where 𝑉is the forward velocity of the vehicle,  𝑉is its derivative, 𝜌ℎ and 𝜌𝑣are the horizontal and vertical curvatures of the designed track respectively, 𝜌𝑡𝑤 is the twist curvature, and 𝜌′ ℎis the spacial derivative of the horizontal curvature with respect to the 𝑠𝑡coordinate. These expressions are of utmost importance for the algorithm proposed in this paper. 3.4. Kinematics of the BF In this work, the TF will refer to the body track frame positioned at the same 𝑋coordinate than the BF’s. Therefore, the relative position of the BF with respect to (w.r.t.) the TF is given by:  𝐫𝑏=⎡ ⎢ ⎢ ⎣ 0 𝑦𝑡,𝑏 𝑧𝑡,𝑏 ⎤ ⎥ ⎥ ⎦ (7) The position and orientation of the BF [𝑂𝑏(𝑋𝑏𝑌𝑏𝑍𝑏)] with respect to the TF is defined with the coordinates: 𝐪𝑡,𝑏 =[𝑠𝑡𝑦𝑡,𝑏 𝑧𝑡,𝑏 𝜑𝑡,𝑏 𝜃𝑡,𝑏 𝜓𝑡,𝑏]𝑇(8) where 𝚽𝑡,𝑏 = [𝜑𝑡,𝑏 𝜃𝑡,𝑏 𝜓𝑡,𝑏]𝑇is a set of three Euler angles that describe the orientation of the BF with respect to the TF. Note the superscript ‘𝑡, 𝑏’ meaning relative position or orientation w.r.t. the track (t) of the body (b). Therefore, Eq. (8) includes only relative coordinates except for the 𝑠𝑡. The position vector 𝐑𝑏of the centre of the BF with respect to the GF, 𝑂𝑏, is given by: 𝐑𝑏=𝐑𝑡+𝐫𝑏=𝐑𝑡+𝐂𝑡 𝐫𝑏(9) where 𝐂𝑡is the rotation matrix that transform the TF components of  𝐫𝑏 into GF components. Once the position of the BF is defined by Eq. (9) and the rotation matrix is obtained from 𝚽𝑡, the calculation of the absolute velocity and acceleration in the TF is given by:   𝐑𝑏=  𝐑𝑡+  𝐫𝑏+ 𝛚𝑡× 𝐫𝑏(10)   𝐑𝑏=  𝐑𝑡+  𝐫𝑏+ 2  𝛚𝑡×  𝐫𝑏+ 𝛂𝑡× 𝐫𝑏+ 𝛚𝑡× (  𝛚𝑡× 𝐫𝑏)(11) where   𝐑𝑡and   𝐫𝑏are the time-derivatives of the position arrays, 𝐑𝑡and 𝐫𝑏, projected to the TF. As TF is a non-fixed frame, it is necessary to account for the frame’s own rotation and add the term  𝛚𝑡× 𝐫𝑏to correct the expression. Cross products of Eqs. (10) and (11) can be expressed with skew-symmetric matrices:   𝐑𝑏=  𝐑𝑡+  𝐫𝑏+  𝛚𝑡 𝐫𝑏(12)   𝐑𝑏=  𝐑𝑡+  𝐫𝑏+ 2   𝛚𝑡  𝐫𝑏+(  𝛂𝑡+  𝛚𝑡  𝛚𝑡) 𝐫𝑏(13) where the matrices   𝛚𝑡and   𝛂𝑡are the skew-symmetric matrices associated with the angular velocity and acceleration of the TF with respect to the GF, respectively. Finally, to express the translational acceleration of the BF in its own local frame, it is necessary to premultiply the Eq. (13) by the rotation matrix (𝐂𝑡,𝑏)𝑇, resulting in:   𝐑𝑏= (𝐂𝑡,𝑏)𝑇  𝐑𝑏= (𝐂𝑡,𝑏)𝑇(  𝐑𝑡+  𝐫𝑏+ 2   𝛚𝑡  𝐫𝑏+(  𝛂𝑡+  𝛚𝑡  𝛚𝑡) 𝐫𝑏)(14) It was important to find this expression because the accelerometer, that is used in this work, measures the absolute acceleration of the body frame in its own local frame. 4. Orientation prediction model No matter what algorithm is used, a prediction model is necessary to estimate the orientation. As stated in Section 1, general purpose applications tend to assume that the signals measured by the IMU are mainly due to gravity while applications of ground vehicles also include translational acceleration in the model. However, none of these approaches make a detailed estimation of the value of that translational acceleration to improve the accuracy of the model. Considering that the IMU measures gravity as well as that absolute translational acceleration in its local frame, it can be stated that: 𝐚𝑖𝑚𝑢 = 𝐠+  𝐑𝑏(15) where 𝐚𝑖𝑚𝑢 stands for the signal measured by the IMU and  𝐠= (𝐂𝑏)𝑇𝐠= (𝐂𝑏)𝑇⎡ ⎢ ⎢ ⎣ 0 0 𝑔⎤ ⎥ ⎥ ⎦ (16) when the IMU is rigidly attached to the BF. This expression can be further expanded using Eq. (14) giving: 𝐚𝑖𝑚𝑢 = (𝐂𝑏)𝑇𝐠+ (𝐂𝑏)𝑇𝐂𝑡(  𝐑𝑡+  𝐫𝑏+ 2   𝛚𝑡  𝐫𝑏+(  𝛂𝑡+  𝛚𝑡  𝛚𝑡) 𝐫𝑏)(17) At this point, the available equations and unknowns can be stated. Considering that the kinematics of the TF can be calculated with Eqs. (3)–(6), the unknowns are the absolute Euler Angles, 𝚽𝑏, that composes the rotation matrix 𝐂𝑏, the position of the body in the TF,  𝐫𝑏, and their time-derivatives. Now, the main hypothesis to make are to assume that 2  𝛚𝑡  𝐫𝑏and   𝐫𝑏are negligible compared with   𝐑𝑡, and that  𝐫𝑏initial value can be approximately measured on-site and its value will not vary much with time. Another very convenient hypothesis is to consider that TF and BF are approximately aligned, therefore 𝐂𝑡,𝑏 = (𝐂𝑡)𝑇𝐂𝑏≈𝐈. These hypotheses are checked in Section 8. This allows to write the Eq. (17) as: 𝐚𝑖𝑚𝑢 ≈ (𝐂𝑏)𝑇𝐠+(  𝐑𝑡+(  𝛂𝑡+  𝛚𝑡  𝛚𝑡) 𝐫𝑏)(18) where a 3 equation - 3 unknown system is reached, with the Euler angles as the sole unknowns when the accelerometer and encoder signals are obtained and  𝐫𝑏magnitude is measured. Once this equation has been setup, different solution methods may be implemented. This equation can be used for the prediction or update model in a Extended Kalman Filter or as function to minimize in a Gradient Descendent Algorithm. These two implementation examples are shown in Section 5. 5. Orientation estimation algorithm In this paper, both an Extended Kalman Filter (EKF) and an optimization algorithm such as the Gradient Descendent Algorithm (GDA) Measurement 184 (2021) 109898 7 E. Briales et al. are proposed for the estimation of 3D orientation. The output of the algorithms is the posing of the BF with respect to the GF, and the inputs are inertial sensor signals and 𝑠𝑡,𝑉, 𝑉and  𝐫𝑏. In this case, the EKF is based on Sabatini’s work [44,45] and the GDA is based on Madgwick’s [26,27]. 5.1. Extended Kalman Filter For the EKF implementation, gyroscope measurements are used for the prediction model and accelerometer measurements are used in the observation model. The prediction model used for the implementation of the EKF is based on the expression that relates the angular velocity with the Euler angles derivatives [44,45]. That is:  𝜔𝑏=𝜓𝑏 𝑘+ 𝜃𝑏 𝑗1+𝜑𝑏 𝑖2(19) Where  𝑘, 𝑗1, 𝑖2are the unit vectors of the global frame and the two intermediate frames (1 and 2, respectively), used in the Euler rotations sequence. If these unit vectors are projected to the BF, it yields:  𝛚𝑏= 𝐆 𝚽𝑏(20) where  𝚽𝑏is the derivative of the Euler angles,  𝛚𝑏is the local angular velocity of the body that it is measured with the gyroscope signals and  𝐆is the matrix whose columns contain the components of the unit vectors in the BF. From now on, the superscript 𝑏will be avoided as the Euler angles are always referred to this body. This equation can be expressed as an ordinary differential equation (ODE) as:  𝚽= 𝐆−1  𝛚=⎡ ⎢ ⎢ ⎣ 1 sin 𝜑tan 𝜃cos 𝜑tan 𝜃 0 cos 𝜑− sin 𝜑 0 sin 𝜑∕ cos 𝜃cos 𝜑∕ cos 𝜃⎤ ⎥ ⎥ ⎦⎡ ⎢ ⎢ ⎣ 𝜔𝑥 𝜔𝑦 𝜔𝑧 ⎤ ⎥ ⎥ ⎦ =𝐟(𝚽)(21) Defining 𝐱(𝑡) = 𝚽(𝑡), the continuous system observation model can be expressed as [53]:  𝐱(𝑡) = 𝐟(𝐱(𝑡), 𝑡) + 𝐰(𝑡)(22) where 𝐰(𝑡)model the continuous-time white noise process with mean zero and intensity 𝐐. This model has to be linearized about a reference state to be implemented. This reference state, 𝐱𝑅, is obtained as solution to the Eq. (22) without the process noise, 𝐰(𝑡)as:  𝐱𝑅(𝑡) = 𝐟(𝐱𝑅(𝑡), 𝑡)(23) An approximate linear model can be presented as:  𝐱(𝑡) −  𝐱𝑅(𝑡) = 𝐟(𝐱(𝑡), 𝑡) + 𝐰(𝑡) − 𝐟(𝐱𝑅(𝑡), 𝑡)(24) The function 𝐟can be linearized about the reference state such as: 𝐟(𝐱(𝑡), 𝑡) ≈ 𝐟(𝐱𝑅(𝑡), 𝑡)+𝜕𝐟(𝐱, 𝑡) 𝜕𝐱||||𝐱=𝐱𝑅 (𝐱−𝐱𝑅)(25) Defining 𝐅matrix as: 𝐅(𝑡) = 𝜕𝐟(𝐱, 𝑡) 𝜕𝐱||||𝐱=𝐱𝑅(𝑡) (26) the Eq. (24) can be linearized such as:  𝐱(𝑡) −  𝐱𝑅(𝑡) = 𝐅(𝑡)(𝐱(𝑡) − 𝐱𝑅(𝑡))+𝐰(𝑡)(27) being Eq. (27) a linear ODE. This system may be split into two different equations that constitutes the prediction step (propagation of the mean and the covariance of the state) and discretized for 𝑘and 𝑘− 1 (new reference state) using the Euler forward integrator such as [54]: 𝐏𝐫𝐞𝐝𝐢𝐜𝐭𝐢𝐨𝐧   𝐱𝑘−1 =𝐟( 𝐱𝑘−1)(28)  𝐱− 𝑘= 𝐱𝑘−1 +  𝐱𝑘−1𝑇𝑠(29)  𝐏− 𝑘=𝐀𝑘−1  𝐏𝑘−1𝐀𝑇 𝑘−1 +∫𝑇𝑠 0(𝑒𝐅𝑘−1(𝑇𝑠−𝜏))𝐐(𝑒𝐅𝑘−1(𝑇𝑠−𝜏))𝑇𝑑𝜏 (30) with the hat () noting that they are estimators of the real state (unknown), the superscript −meaning a priori estimations, 𝑇𝑠being the sampling time, with 𝐀𝑘−1 defined as: 𝐀𝑘−1 =𝑒𝐅𝑘−1𝑇𝑠(31) 𝐅𝑘−1 being the matrix resulting of the linearization of the process model such as: 𝐅𝑘−1 =𝜕𝐟 𝜕𝐱||||𝐱=𝐱𝑘−1 (32) and 𝐐𝑘−1 is the process noise covariance matrix supposed proportional to the sampling period, 𝑇𝑠, such as 𝐐𝑘−1 =𝑘𝑇𝑠𝐈, with 𝐈the 3×3 identity matrix and 𝑘the proportionality constant. The process noise covariance matrix is propagated with the integral in Eq. (30) that can be computed using Van Loan’s Method [54,55]. However, the update model uses the kinematic-specific equations presented in Eq. (17). The observation model can be expressed as: 𝐳𝑘=𝐚𝑖𝑚𝑢 𝑘−(  𝐑𝑡 𝑘+ (  𝛂𝑘+  𝛚𝑘  𝛚𝑘) 𝐫𝑏 𝑘)= (𝐂𝑏 𝑘)𝑇𝐠=𝐡(𝐱𝑘)(33) This equation has the form of 𝐳𝑘=𝐡(𝐱𝑘).𝐳𝑘gathers the signals measured by the accelerometer without the estimated translational acceleration at each instant of time and 𝐡is the measurement function. Therefore, the update equations can implemented be as: 𝐔𝐩𝐝𝐚𝐭𝐞 𝐯𝑘=𝐳𝑘−𝐡( 𝐱− 𝑘)(34) 𝐒𝑘=𝐇𝑘 𝐏− 𝑘𝐇𝑇 𝑘+𝐑𝑘(35) 𝐊𝑘= 𝐏− 𝑘𝐇𝑇 𝑘𝐒−1 𝑘(36)  𝐱𝑘= 𝐱− 𝑘+𝐊𝑘𝐯𝑘(37)  𝐏𝑘= 𝐏− 𝑘−𝐊𝑘𝐒𝑘𝐊𝑇 𝑘(38) where 𝐯𝑘, known as the innovation, is the absolute error between real measurements and the values estimated by the measurement function; 𝐑𝑘is the sensor noise covariance matrix; 𝐒𝑘gathers the uncertainty of prediction and measurement; 𝐊𝑘is the so-called Kalman gain and 𝐇𝑘 is the linearization of Eq. (33) such as: 𝐇𝑘=𝜕𝐡 𝜕𝐱||||𝐱=𝐱− 𝑘 (39) It can be observed that this model is similar to that used for ground vehicles by Klier [39] and others. Nevertheless, in this paper, the expression for the translational acceleration is presented in much more detail, considering the design geometry of the guideway is known in advance. The result of this implementation is shown in Section 6. 5.2. Gradient descendent algorithm As suggested by Madgwick et al. [26,27], this method provides the optimal orientation of the sensor frame such that the measurement of the accelerometer coincides with the gravity vector 𝐠. The function to minimize in that paper is: 𝐟𝐌𝐚𝐝𝐠𝐰𝐢𝐜𝐤(𝚽𝑏,𝐚𝑖𝑚𝑢) = 𝐚𝑖𝑚𝑢 − 𝐠=𝐚𝑖𝑚𝑢 − (𝐂𝑏)𝑇𝐠(40) Naturally, the accelerometer also measures the acceleration of the sensor. Therefore, a more accurate function to minimize includes the translational acceleration term as given in Section 4, as follows: 𝐟(𝚽𝑏,𝐚𝑖𝑚𝑢) = 𝐚𝑖𝑚𝑢 − (𝐂𝑏)𝑇𝐠−  𝐑𝑏(41) where   𝐑𝑏is the absolute translational acceleration of the BF with respect to the GF projected to the BF. It must be noted that these translational accelerations are, in general applications, unknown. Consequently, the common approach is to disregard this term. However, in railway applications this acceleration can be well approximated if the track geometry and the forward velocity of the vehicle are known, Measurement 184 (2021) 109898 8 E. Briales et al. Fig. 5. Experimental setup of the IMU sensor in the railway axle (body frame). as it was previously stated. Considering Eq. (13) and the hypotheses made in previous sections, the minimization function of Eq. (41) can be rewritten as: 𝐟(𝚽𝑏,𝐚𝑖𝑚𝑢) = 𝐚𝑖𝑚𝑢 − (𝐂𝑏)𝑇𝐠−(  𝐑𝑡+(  𝛂𝑡+  𝛚𝑡  𝛚𝑡) 𝐫𝑏)(42) Once the minimization function is set, the GDA algorithm can be used with the implementation detailed in Madgwick’s works [26, 27]. The results of this algorithm are shown in Section 6. More recently, Madgwick et al. [24] published a new orientation estimation article showing promising results. However, that article is mainly focused on improving performance under magnetic distortions while the method proposed in this paper does not use magnetometer data. 6. Experimental setup and rides The experiments were made in a light-metro train property of Metro de Sevilla riding on its operational track at different speeds. This lightmetro vehicle uses independently rotating wheels (IRW). The use of IRW allows the installation of the inertial sensor on the wheelset axle (see Fig. 5). Therefore, this solid is considered the body frame whose origin is coincident with the sensor frame’s. The complete track geometry measurement system (TGMS) also included a computer vision system to measure the displacements of the railheads with respect to the body frame. However, the output data of these sensors is not used for this orientation estimation algorithm and therefore the camera system shown in Fig. 5 is not described. The IMU used was a AHRS 3DM-GX4®-25 of LORD MicroStrain that provided the measurements of acceleration, angular velocity, magnetometer and Euler angles at 250 Hz. However, the magnetometer signal was finally not considered for the proposed orientation estimation algorithms because this sensor is most likely subjected to strong magnetic distortions due to the industrial environment and the underground rides. In Table 1, the main technical specifications of the sensor are shown. The data of this IMU was collected with a ad-hoc data acquisition system based on a microcontroller which gathered the measurements with a RS232 protocol and send them to a computer with a serial communication protocol (see Fig. 6). The different parameter values used for the implementation of the algorithm where the following: The measured value of  𝐫𝑏was [0,0,0.6]𝑇metres for the initial position. For the Madgwick’s algorithm implementation, the value of the estimated mean zero gyroscope Table 1 Technical specifications of sensor AHRS 3DM-GX4®-25. Accelerometer Gyroscope Magnetometer Measurement range ±5 g ±300◦/sec ±2.5 Gauss Non-linearity error ±0.03% fsa±0.03% fs ±0.4% fs Resolution <0.1 mg 0.008◦/sec – Bias instability ±0.04 mg 10◦/hr – Initial bias error ±0.002 g ±0.05◦/sec ±0.003 Gauss Noise density 100 μg/√Hz 0.005◦/sec/√Hz 100 μGauss/√Hz Alignment error ±0.05◦±0.05◦±0.05◦ Acquisition frequencyb4 kHz 4 kHz 50 Hz max afull scale, % of error with respect to the real absolute value. bThis acquisition frequency is much higher than the sampling frequency to average the output signal and reduce the noise. Fig. 6. Data acquisition system hardware used during the experiments. measurement error, 𝜔𝐸𝑟𝑟𝑜𝑟, was 0.5◦∕𝑠. For the EKF algorithm implementation, the values of the covariance matrices were initially set to 𝐐𝑘−1 =𝑘𝑇𝑠𝐈= 1 ⋅10−2𝐈with 𝑘= 2.5and 𝑇𝑠= 4 ms, and 𝐑𝑘=𝜎2 𝑎𝑐𝑐𝐈 with 𝜎𝑎𝑐𝑐 = 1.5⋅10−2 m∕s2. However, they were adjusted to new values to obtain better results and the final settings implemented were 𝐐𝑘−1 = 1 ⋅10−3𝐈and 𝐑𝑘= (1.5)2⋅102𝐈. This decreases the uncertainty of the model and increases the uncertainty of the measurement, so the Measurement 184 (2021) 109898 9 E. Briales et al. Fig. 7. Real instantaneous velocities of the different rides. Table 2 Comparative between the characteristics of the different rides in which measurements were taken. Measurements 01/03 02/04 05 06 Identification go20A/go20B back20A/back20B go70 back70 Origin station Condequinto Olivar Condequinto Olivar Origin PK 15+900 18+00 15+900 18+000 Destination PK 18+000 15+900 18+000 15+900 Destination station Olivar Condequinto Olivar Condequinto Length (m) 2100 −2100 2100 −2100 Ride direction forwards backwards forwards backwards Velocity low low commercial commercial Velocity setpoint (km/h) 20 −20 70 −70 Mean velocity (km/h) 21/18 −21/−21 42 −31.5 algorithm mainly uses the prediction for the orientation estimation but the drift is still corrected by the update step. In order to test the validity of the algorithm proposed in this paper, different rides where made to collect inertial and forward velocity data. The different riding tests are summarized in Table 2. The forward velocity has an important influence in the estimation of the orientation because, along with the geometry, velocity has direct influence on accelerations and angular velocities of the track frame, as it can be seen in Eqs. (3)–(6). Therefore, low-speed and high-speed rides will be considered independently to properly understand the differences. The instantaneous velocities of these rides are shown in Fig. 7, where the velocity in the backward direction is considered as negative. The vertical and horizontal projection of the track centreline covered during the experiments are shown in Figs. 8 and 9, and the curvatures of the sections are shown in Fig. 10. This route contains smooth curves, sharp curves and straight sections, what permits to evaluate the behaviour of the proposed algorithms for different track section geometries. The sensor signal and estimated track frame kinematics must be referenced to the same point in the track. For this reason, position of the vehicle along the track at each instant is necessary. The authors have developed their own odometry algorithm that corrects the position of the railway by comparing the vertical angular velocity measured when the vehicle is entering or exiting a curve. This method allows to correct the errors in odometry reaching an accuracy in the order of magnitude of 1 metre. This algorithm is outside the scope of this paper but the readers can refer to Escalona’s work for further explanation [47]. The results of proposed orientation estimation method for the rides given in Table 2 are shown in Section 7. 7. Comparative The proposed estimation methods have been implemented using the presented Extended Kalman Filter and the corrected Madgwick’s algorithm. In this section, results have been compared to the orientation given by commercial IMU for different rides. 7.1. EKF and Madgwick comparative for Test 01 As it can be seen in Fig. 11 for the ride go20A, the improvements of the orientation estimation when introducing the kinematics of the track frame are noticeable. Not only the yaw curve has a reasonable shape, but the roll and pitch do follow the reference value of these angles properly. When comparing EKF implementation with Madgwick’s, it is also worth noting that the results are very similar for both algorithms. As gravity can only be used to calculate roll and pitch (heading does not modify local gravity in IMU axes) the calculation of yaw angle is exactly the same for EKF and GDA and therefore is subjected to the same drift issue. This makes the EKF and Madgwick’s values for this angle appear overlapped in the figures. At this point, the different behaviour of the Euler angles obtained with GDA and EKF are clear. In the case of the roll angle, while the GDA properly follow the track reference angles, the EKF fails to obtain the proper value both in curves and in straights where the angle does not return to 0. In the case of the pitch angle, both methods are very similar, having the same overall error (around ±1◦). The difference between both methods for roll angles requires a deeper reflection. While the EKF is usually the algorithm that grants best results, that was not the case for the authors. At first, it was observed that Euler angles results had little response to big modifications the EKF parameters (𝐐𝑘−1,𝐑𝑘). Which means that the sensitivity of the algorithm to these parameters was very low. To improve the results, the relative weight between the Measurement 184 (2021) 109898 16 E. Briales et al. Fig. 21. Rides go20A +back20A +go20B +back20B +go70 +back70. Madgwick orientation estimations. Fig. 22. Rides go20A +back20A +go20B +back20B +go70 +back70. EKF orientation estimations. several parameters to tune-up. This is considered to be mainly due to the complication of establishing proper Kalman uncertainty parameters, while in Madgwick’s algorithm optimal parameters values were already provided by the author. Once the prediction model is more realistic, the accuracy of the orientation estimation may be increased improving the different algorithm implementations. Moreover, the use of certain constrains based on known railway kinematics and track kinematics may also be included in the different algorithms. Finally, the use of certain reliability tests to modify the weight of the different measurements when introduced into the algorithm would also improve the final results of the algorithm. CRediT authorship contribution statement Eduardo Briales: Conceptualization, Software, Investigation, Data curation, Writing – original draft, Visualization. Pedro Urda: Investigation. José L. Escalona: Conceptualization, Methodology, Investigation, Supervision, Writing – review and editing, Validation. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The first author thanks the support of the Spanish Ministry of Economy, Industry and Competitiveness under the project reference DI-15-07658. This support is gratefully acknowledged. The second and third authors thank the support of the Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía, through the program ‘Programa Operativo FEDER 2014- 2020’, awarded to the University of Seville, financed with the European Regional Development Fund (FEDER), under Project US-1257665. This work could have not been possible without the close collaboration with Metro de Sevilla and its Rolling Stock and Infrastructure Manager, Enrique Martínez Asensio. Measurement 184 (2021) 109898 17 E. Briales et al. References [1] H. Fourati, N. Manamanni, L. Afilal, Y. Handrich, Sensors-based data fusion solution design for 3D motion estimation with application in Bio-logging, Int. J. Sci. Tech. Autom. Control Comput. Eng. 3 (2) (2009) 1012–1031. [2] L. Hu, W. Yang, J. He, H. Zhou, Z. Zhang, X. Luo, R. Zhao, L. Tang, P. Du, Roll angle estimation using low cost MEMS sensors for paddy field machine, Comput. Electron. Agric. 158 (February) (2019) 183–188, http://dx.doi.org/10. 1016/j.compag.2019.01.010. [3] P. Huang, Z. Zhang, X. Luo, J. Zhang, P. Huang, Attitude estimation of agricultural implements based on quaternion and complementary filter, IFAC- PapersOnLine 51 (17) (2018) 837–842, http://dx.doi.org/10.1016/j.ifacol.2018. 08.090. [4] C. Miao, H. Chu, J. Cao, Z. Sun, R. Yi, Steering angle adaptive estimation system based on GNSS and MEMS gyro, Comput. Electron. Agric. 153 (July) (2018) 196–201, http://dx.doi.org/10.1016/j.compag.2018.08.019. [5] S.P. Driessen, N.H. Janssen, L. Wang, J.L. Palmer, H. Nijmeijer, Experimentally validated extended Kalman filter for UAV state estimation using low-cost sensors, IFAC-PapersOnLine 51 (15) (2018) 43–48, http://dx.doi.org/10.1016/j.ifacol. 2018.09.088. [6] W. Alqaisi, J. Ghommam, A. Alazzam, M. Saad, V. Nerguizian, Three-loop uncertainties compensator and sliding mode quadrotor control, Comput. Electr. Eng. 81 (2020) 106507, http://dx.doi.org/10.1016/j.compeleceng.2019.106507. [7] T.H. Bryne, R.H. Rogne, T.I. Fossen, T.A. Johansen, Attitude and heave estimation for ships using MEMS-based inertial measurements, IFAC-PapersOnLine 49 (23) (2016) 568–575, http://dx.doi.org/10.1016/j.ifacol.2016.10.496. [8] M. Aldimirov, R. Arnaudov, Method for automated reconstruction of a car’s path during crash from GPS/INS data using a Kalman filter, Adv. Eng. Softw. 115 (November 2017) (2018) 386–390, http://dx.doi.org/10.1016/j.advengsoft.2017. 10.009. [9] W.J. Scholte, V. Rodrigo Marco, H. Nijmeijer, Experimental validation of vehicle velocity, attitude and IMU bias estimation, IFAC-PapersOnLine 52 (8) (2019) 118–123, http://dx.doi.org/10.1016/j.ifacol.2019.08.058. [10] J. Vaganay, M.J. Aldon, Attitude estimation for a vehicle using inertial sensors, Control Eng. Pract. 2 (2) (1994) 281–287, http://dx.doi.org/10.1016/0967- 0661(94)90209-7. [11] X. Xu, X. Tian, L. Zhou, Y. Li, A decision-tree based multiple-model UKF for attitude estimation using low-cost MEMS marg sensor arrays, Meas. J. Int. Meas. Confed. 135 (2019) 355–367, http://dx.doi.org/10.1016/j.measurement.2018.11. 062. [12] M. Anderle, S. Čelikovský, Sensor fusion for simple walking robot using low-level implementation of extended Kalman filter, IFAC-PapersOnLine 51 (13) (2018) 43–48, http://dx.doi.org/10.1016/j.ifacol.2018.07.252. [13] X. Li, W. Chen, C. Chan, B. Li, X. Song, Multi-sensor fusion methodology for enhanced land vehicle positioning, Inf. Fusion 46 (March 2017) (2019) 51–62, http://dx.doi.org/10.1016/j.inffus.2018.04.006. [14] T. Konrad, J.J. Gehrt, J. Lin, R. Zweigel, D. Abel, Advanced state estimation for navigation of automated vehicles, Annu. Rev. Control 46 (2018) 181–195, http://dx.doi.org/10.1016/j.arcontrol.2018.09.002. [15] G.A. Aydemir, A. Saranli, Characterization and calibration of MEMS inertial sensors for state and parameter estimation applications, Meas. J. Int. Meas. Confed. 45 (5) (2012) 1210–1225, http://dx.doi.org/10.1016/j.measurement. 2012.01.015. [16] V.R. Marco, J. Kalkkuhl, J. Raisch, EKF for simultaneous vehicle motion estimation and IMU bias calibration with observability-based adaptation, in: Proceedings of the American Control Conference, Vol. 2018, 2018, pp. 6309–6316, doi:10.23919/ACC.2018.8431245. [17] D. Tcherniak, M. Schwaab, On a method for finding position and orientation of accelerometers from their signals, Mech. Syst. Signal Process. 140 (2020) 106662, http://dx.doi.org/10.1016/j.ymssp.2020.106662. [18] I. Ahmad, A. El Hadri, L. Benziane, A. Benallegue, Globally asymptotic attitude estimation for accelerated aerial vehicles, Aerosp. Sci. Technol. 84 (2019) 1175–1181, http://dx.doi.org/10.1016/j.ast.2018.06.008. [19] M.D. Hua, Attitude estimation for accelerated vehicles using GPS/INS measurements, Control Eng. Pract. 18 (7) (2010) 723–732, http://dx.doi.org/10.1016/j. conengprac.2010.01.016. [20] D.A. Aligia, B.A. Roccia, C.H. De Angelo, G.A. Magallán, G.N. González, An orientation estimation strategy for low cost IMU using a nonlinear luenberger observer, Meas. J. Int. Meas. Confed. (April) (2020) 108664, http://dx.doi.org/ 10.1016/j.measurement.2020.108664. [21] B. Kada, K. Munawar, M.S. Shaikh, M.A. Hussaini, U.M. Al-Saggaf, UAV Attitude estimation using nonlinear filtering and low-cost mems sensors, IFAC- PapersOnLine 49 (21) (2016) 521–528, http://dx.doi.org/10.1016/j.ifacol.2016. 10.655. [22] H. Benzerrouk, A. Nebylov, H. Salhi, Quadrotor UAV state estimation based on high-degree Cubature Kalman filter, IFAC-PapersOnLine 49 (17) (2016) 349–354, http://dx.doi.org/10.1016/j.ifacol.2016.09.060. [23] M. Kok, T.B. Schön, A fast and robust algorithm for orientation estimation using inertial sensors, IEEE Signal Process. Lett. 26 (11) (2019) 1673–1677, http://dx.doi.org/10.1109/LSP.2019.2943995,arXiv:1910.00463. [24] S.O. Madgwick, S. Wilson, R. Turk, J. Burridge, C. Kapatos, R. Vaidyanathan, An extended complementary filter for full-body MARG orientation estimation, IEEE/ASME Trans. Mechatronics 25 (4) (2020) 2054–2064, http://dx.doi.org/ 10.1109/TMECH.2020.2992296. [25] M. Abolhasani, M. Rahmani, Robust deterministic least-squares filtering for uncertain time-varying nonlinear systems with unknown inputs, Systems Control Lett. 122 (2018) 1–11, http://dx.doi.org/10.1016/j.sysconle.2018.09.005. [26] S.O. Madgwick, An efficient orientation filter for inertial and inertial/magnetic sensor arrays, Report X-Io and University of ..., 2010, p. 32, http://dx.doi.org/ 10.1109/ICORR.2011.5975346. [27] S.O. Madgwick, A.J. Harrison, R. Vaidyanathan, Estimation of IMU and MARG orientation using a gradient descent algorithm, in: IEEE International Conference on Rehabilitation Robotics, 2011, pp. 1–7, doi:10.1109/ICORR.2011.5975346. [28] S. Wilson, H. Eberle, Y. Hayashi, S.O. Madgwick, A. McGregor, X. Jing, R. Vaidyanathan, Formulation of a new gradient descent MARG orientation algorithm: Case study on robot teleoperation, Mech. Syst. Signal Process. 130 (2019) 183–200, http://dx.doi.org/10.1016/j.ymssp.2019.04.064. [29] S. Mansoor, U.I. Bhatti, A.I. Bhatti, S.M.D. Ali, Improved attitude determination by compensation of gyroscopic drift by use of accelerometers and magnetometers, Meas. J. Int. Meas. Confed. 131 (2019) 582–589, http://dx.doi.org/10.1016/j. measurement.2018.08.067. [30] D. Li, J. Zhou, Y. Liu, Recurrent-neural-network-based unscented Kalman filter for estimating and compensating the random drift of MEMS gyroscopes in real time, Mech. Syst. Signal Process. 147 (2021) 107057, http://dx.doi.org/10.1016/ j.ymssp.2020.107057. [31] H. Nourmohammadi, J. Keighobadi, Fuzzy adaptive integration scheme for lowcost SINS/GPS navigation system, Mech. Syst. Signal Process. 99 (2018) 434–449, http://dx.doi.org/10.1016/j.ymssp.2017.06.030. [32] X. Xu, Y. Sun, X. Tian, L. Zhou, Y. Li, A proposed attitude estimator with reliability test criteria for sensor data fusion, Meas. J. Int. Meas. Confed. 150 (2020) 107046, http://dx.doi.org/10.1016/j.measurement.2019.107046. [33] J.G. Rueda-Escobedo, J.A. Moreno, Delayed Kalman-bucy observer for a class of LTV systems with delayed measurements, IFAC-PapersOnLine 51 (25) (2018) 12–17, http://dx.doi.org/10.1016/j.ifacol.2018.11.074. [34] A. Katriniok, D. Abel, Adaptive EKF-based vehicle state estimation with online assessment of local observability, IEEE Trans. Control Syst. Technol. 24 (4) (2016) 1368–1381, http://dx.doi.org/10.1109/TCST.2015.2488597. [35] F. Zhu, W. Zhou, Y. Zhang, R. Duan, X. Lv, X. Zhang, Attitude variometric approach using DGNSS/INS integration to detect deformation in railway track irregularity measuring, J. Geod. 93 (9) (2019) 1571–1587, http://dx.doi.org/10. 1007/s00190-019-01270-w. [36] M. Li, H. Yu, X. Zheng, A.I. Mourikis, High-fidelity sensor modeling and self-calibration in vision-aided inertial navigation, in: Proceedings - IEEE International Conference on Robotics and Automation, 2014, pp. 409–416, doi: 10.1109/ICRA.2014.6906889. [37] E. Zappa, R. Liu, L. Trainelli, A. Ferrario, P. Cordisco, M. Terraneo, R. Grassetti, M. Redaelli, Laser and vision-based measurements of helicopter blade angles, Meas. J. Int. Meas. Confed. 118 (August 2017) (2018) 29–42, http://dx.doi.org/ 10.1016/j.measurement.2017.12.037. [38] L. Zhu, Z. Yu, Instantaneous position and pose measurements of moving vehicles with applications to railway infrastructure monitoring, in: 2009 IEEE Intrumentation and Measurement Technology Conference, I2MTC 2009, 2009, pp. 553–556, doi:10.1109/IMTC.2009.5168511. [39] W. Klier, A. Reim, D. Stapel, Robust estimation of vehicle sideslip angle-an approach w/o vehicle and tire models, SAE Tech. Pap. 2008 (724) (2008) http://dx.doi.org/10.4271/2008-01-0582. [40] G. Dissanayake, S. Sukkarieh, E. Nebot, H. Durrant-Whyte, The aiding of a lowcost strapdown inertial measurement unit using vehicle model constraints for land vehicle applications, IEEE Trans. Robot. Autom. 17 (5) (2001) 731–747, http://dx.doi.org/10.1109/70.964672. [41] H.F. Grip, L. Imsland, T.A. Johansen, J.C. Kalkkuhl, A. Suissa, Estimation of road inclination and bank angle in automotive vehicles, in: Proceedings of the American Control Conference, 2009, pp. 426–432, doi:10.1109/ACC.2009. 5159912. [42] Q. Yuan, E. Asadi, Q. Lu, G. Yang, I.M. Chen, Uncertainty-based IMU orientation tracking algorithm for dynamic motions, IEEE/ASME Trans. Mechatronics 24 (2) (2019) 872–882, http://dx.doi.org/10.1109/TMECH.2019.2892069. [43] A. Makni, H. Fourati, A.Y. Kibangou, Adaptive Kalman filter for MEMS-IMU based attitude estimation under external acceleration and parsimonious use of gyroscopes, in: 2014 European Control Conference, ECC 2014, 2014. pp. 1379–1384, doi:10.1109/ECC.2014.6862535. [44] A. Sabatini, Quaternion-based extended Kalman filter for determining orientation by inertial and magnetic sensing, IEEE Trans. Biomed. Eng. 53 (7) (2006) 1346– 1356, http://dx.doi.org/10.1109/TBME.2006.875664,http://ieeexplore.ieee.org/ document/1643403/. Measurement 184 (2021) 109898 18 E. Briales et al. [45] A.M. Sabatini, Estimating three-dimensional orientation of human body parts by inertial/magnetic sensing, Sensors 11 (2) (2011) 1489–1525, http://dx.doi.org/ 10.3390/s110201489. [46] J.F. Aceituno, R. Chamorro, S. Muñoz, J.L. Escalona, An alternative procedure to measure railroad track irregularities. application to a scaled track, Meas. J. Int. Meas. Confed. 137 (2019) 417–427, http://dx.doi.org/10.1016/j.measurement. 2019.01.025. [47] J.L. Escalona, P. Urda, S. Muñoz, A track geometry measuring system based on multibody kinematics, inertial sensors and computer vision, Sensors 21 (3) (2021) 683, http://dx.doi.org/10.3390/s21030683,https://www.mdpi.com/ 1424-8220/21/3/683. [48] T.H. Riehle, S.M. Anderson, P.A. Lichter, N.A. Giudice, S.I. Sheikh, R.J. Knuesel, D.T. Kollmann, D.S. Hedin, Indoor magnetic navigation for the blind, in: Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS, IEEE, 2012, pp. 1972–1975, http: //dx.doi.org/10.1109/EMBC.2012.6346342. [49] N. Yadav, C. Bleakley, Accurate orientation estimation using AHRS under conditions of magnetic distortion, Sensors (Switzerland) 14 (11) (2014) 20008–20024, http://dx.doi.org/10.3390/s141120008. [50] E. Palermo, S. Rossi, F. Patanè, P. Cappa, Experimental evaluation of indoor magnetic distortion effects on gait analysis performed with wearable inertial sensors, Physiol. Meas. 35 (3) (2014) 399–415, http://dx.doi.org/10.1088/0967- 3334/35/3/399. [51] J.L. Escalona, J.F. Aceituno, Multibody simulation of railway vehicles with contact lookup tables, Int. J. Mech. Sci. 155 (January 2018) (2019) 571–582, http://dx.doi.org/10.1016/j.ijmecsci.2018.01.020. [52] S. Muñoz, J.F. Aceituno, P. Urda, J.L. Escalona, Multibody model of railway vehicles with weakly coupled vertical and lateral dynamics, Mech. Syst. Signal Process. 115 (2019) 570–592, http://dx.doi.org/10.1016/j.ymssp.2018.06.019. [53] D. Morrell, Extended Kalman filter lecture notes, in: EEE 581 - Spring 1997, 1997. [54] J. Hartikainen, S. Särkkä, Optimal filtering with Kalman filters and smoothers – a manual for Matlab toolbox EKF / UKF, J. Interprofessional Care 25 (December 2013) (2007) 1–57. [55] C. Van Loan, Computing integrals involving the matrix exponential, IEEE Trans. Automat. Control 23 (3) (1978) 395–404, http://dx.doi.org/10.1109/TAC.1978. 1101743,http://ieeexplore.ieee.org/document/1101743/.