Sensor Saturation Compensated Smoothing Algorithm for Inertial Sensor Based Motion Tracking

In this paper, a smoothing algorithm for compensating inertial sensor saturation is proposed. The sensor saturation happens when a sensor measures a value that is larger than its dynamic range. This can lead to a considerable accumulated error. To compensate the lost information in saturated sensor data, we propose a smoothing algorithm in which the saturation compensation is formulated as an optimization problem. Based on a standard smoothing algorithm with zero velocity intervals, two saturation estimation methods were proposed. Simulation and experiments prove that the proposed methods are effective in compensating the sensor saturation.


Introduction
In motion tracking, there are many ways to estimate the trajectory of a moving object. Moving objects can be tracked accurately by using visual devices such as camera systems [1,2]. In [1], with 40 landmarks attached on the body, Hong et al. used the Eagle Digital motion capture system with seven cameras to analyze 14 angles and one ratio of gait features. In [2], Lee and Grimson investigated person identification and gender classification based on moments computed from the silhouette of walking people. However, these camera systems are limited in their setup ranges and sometimes have high implementation costs. Moreover, the angle views of the cameras are also limited and they are OPEN ACCESS easily affected by illumination. Due to these reasons, for long distance or outdoor measurements, motion tracking based on camera systems seems to be a difficult task.
To avoid the mentioned disadvantages, inertial measurement units (IMU) can be used instead as wearable devices. IMUs are widely used due to their small size and low cost. With the development of the technology, IMUs are now becoming more accurate. In [3], Tadano et al. proposed a method using quaternion calculations from seven sensor units consisting of a tri-axial acceleration and gyro sensors. The quaternions, which are computed from the sensors attached on limbs and waist, are used in a gait wire frame model to generate the gait animation. To increase the accuracy, the IMUs are used with other aiding devices such as cameras in [4,5] or force sensors in [6].
However, IMUs still have their own limitations such as susceptibility to noise and limited dynamic range. The accuracy of inertial sensor-based estimation can be improved by using zero velocity updates as in [7], or taking advantage of the relative position and attitude of multiple sensors [8,9]. In [7], a robot arm control with an automatic calibration function based on inertial sensors is proposed. The authors state that the drift of the sensors is clearly removed by applying zero velocity updates. In [8], Helten et al. introduce another way which takes into account the relative position and attitude of multiple sensors to improve the accuracy of human motion estimation. Using the same method, Tao et al. [9] estimate limb movement. They also use the limb biomechanical model characteristics to provide constraints for sensors' relations. However, the other drawback of the inertial sensor, the saturation, has not been considered. Saturation is a state in which the signal that needs to be measured is larger than the dynamic range of the sensor. When that happens, the output of the sensor becomes the limiting value of the sensor range. This induces a considerable error between the true and estimated values during motion tracking. In this paper, we propose two methods for estimating the sensor saturation. Assuming that the motion is between not moving intervals (called zero velocity intervals), we formulate a saturation estimation problem as an optimization problem by modifying a smoother used in [10,11]. In the proposed methods, some state constraints mentioned in [12] could be added to increase the accuracy of the algorithm.
The paper is organized in five main sections and a conclusion. Section 2 points out the problem formulation. In Section 3, a standard smoothing algorithm with zero velocity intervals is described in detail. Sections 4 and 5 propose some methods for sensor saturation estimation. Some experiments to verify the proposed methods are given in Section 6. The last section concludes the paper.

Problem Formulation
We consider a moving object case where an IMU is attached on the object. There are two coordinate systems in this paper: the navigation coordinate frame and the body coordinate frame. The z axis of the navigation coordinate frame coincides with the local vertical. The choice of x axis is arbitrary. The body coordinate frame is defined as a frame with three axes coinciding with the three axes of the inertial measurement unit. The subscripts b (body) and n (navigation) are used to emphasize that a vector or matrix belongs to the body or navigation coordinate frame, respectively.
Our goal is to estimate attitude (expressed using the quaternion), position and velocity of the object from the sensor data. Let q q q q q R  be a quaternion representing the rotation relationship between the navigation coordinate frame and the body coordinate frame. Let   33 C q R   be the rotation matrix corresponding to the quaternion q [13], and 3 rR  and 3 vR  be the position and velocity of the object, respectively. We have the following basic equations [14]: where 3 n aR  and 3 b aR  are the acceleration made by forces other than gravitational field in the navigation coordinate frame and body coordinate frame, respectively. The symbol  is defined by: The inertial measurement unit used in this paper consists of three axis gyroscopes and accelerometers. Let 3 g yR  be the gyroscope output and 3 a yR  be the accelerometer output. They satisfy the following relationship [15]: yn y a C q g n where g n and a n are zero mean white Gaussian sensor noises with covariances   The symbol g denotes the gravitational acceleration. It is assumed the sensor bias is already compensated using the standard calibration algorithm [16].
In summary, our goal is to estimate the quaternion, velocity and position of a moving object using the accelerometer and gyroscope data where there could be sensor saturation. Since a smoother is used instead of a filter, we note that the proposed method is offline analysis of attitude and position.

Standard Smoothing Algorithm with Zero Velocity Intervals
In this section, a standard smoothing algorithm is formulated in the quadratic optimization problem. A general method of formulating the smoothing problem in the optimization problem is given in [10]. In this section, we apply the result in [10] to an attitude and position estimation problem, where there are zero velocity intervals. Sensor saturation compensation is not considered in this section and will be discussed in Section 4.
We assume that the motion is a short movement consisting of a moving interval and two not moving intervals (see Figure 1). This type of movement can be found in walking [4,5], golf swing [17] and so on. The sensor data are sampled with the sampling period T . It is assumed that there are total N sampling sensor data and the moving interval starts at 1 kT and stop at 2 kT( 12 k k N  ). The rest are the zero velocity intervals. We use the subscript k to describe a variable is expressed in the discrete time k . For example, a gyroscope output data at time k is denoted by , gk y .

Figure 1.
Standard inertial navigation algorithm with zero velocity correction.
Using the gyroscope output g y in Equation (1), the q , v and r can be estimated by the following equation ("  " denotes for estimation): The initial values of position 0 r and velocity 0 v are assumed to be zero due to the fact that the object is not moving at the not moving period and the origin of the navigation coordinate frame coincides with the starting point of the object. The initial quaternion 0 q is obtained from the accelerometer data using the following (note that 0 b a  during the zero velocity intervals): The heading in 0 q is not determined and can be chosen arbitrarily. Since the noise terms are included in g y and a y in Equation (3), q , v and r are different from the true q , v and r . The error in q , v and r are estimated using a smoothing algorithm. Thus a smoothing algorithm is not used for directly estimating q , v and r but for estimating errors in q , v and r . Once the errors are estimated, we can update q , v and r to obtain more accurate estimates. A multiplicative error is used for the error in q [18]. A small error in q is denoted by e q . e q is assumed to satisfy the following: ˆe q q q  (4) or in matrix expression: The estimation error in ˆ, qv and r can be expressed using the following state: This   xt is estimated using a smoother. To do that, we derive a differential equation for   xt.
From the assumption that e q and g y   is small, we obtain the follow (see [19] for the derivation): Since the sensor sampling period is T, Equation (8) is discretized with the sampling period T [20]: Since no external sensors other than the inertial sensor are used, there is no physical measurement during the motion. However, the fact that the velocity is zero during the not moving period can be used as a virtual measurement. In motion analysis, it is assumed that the object is not moving when the gyroscope and the variation of the accelerometer are smaller than threshold values for some specified time. Thus there will be a chance that a moving interval is detected as a zero velocity interval. To reflect this fact, a small noise , vk n is added in the following equation: The noise , vk n at time k is modeled as a Gaussian white noise with the covariance . Equation 10 can be rewritten in form of following expression in discrete time at where   . Equation (10) can be used during the not moving interval.
We introduce a set m Z which consist of the discrete time indices belonging to the zero velocity intervals. That is, if m kZ  , then we can use the Equation (11). A smoother problem to estimate k x can be formulated as the following optimization problem [10,11]: Find k x and k w for 0 kN  that minimize: It is assumed that the initial value 0 x and the initial covariance error 0 P are given. The method to choose these initial values will be discussed Section 5. The Equation (12)  x A x w   into Equation (12), we can remove the variable k w from the optimization problem: Let the optimization variable x be defined by , then the matrix form of the optimization will be: where (13). Note that Equation (14) is a quadratic function of x , which can be computed efficiently using the quadratic optimization method [21]. Minimizing Equation (14) will provide a set of estimation error. From these values, ˆ, qv andr can be updated using Equations (4) and (7). Let the minimum solution to the problem. Equation (14) In Equation (14), constraints can be easily added to improve the accuracy. For example, in gait analysis, while walking on a plane that assumed to be parallel with the xy plane of the local navigation coordinate frame, we can make use of zero z axis position as following:

Sensor Saturation Estimation
The sensor saturation occurs when the measured values are over the sensors' dynamic ranges ( Figure 2). While it is difficult to know sensor noise values ( g n and a n in Equation (7)), it is not difficult to know when the sensor saturation occurs.
Each sensor has its own range of measurement. When the measured values are larger than the measurement limit, the saturation happens. Even if the saturation happens in a short time, it could lead to a large accumulated error since the lost information is in a large magnitude data area. The data loss due to the saturation has a considerable influence to the result, especially when the integration is used in data processing.
The saturation can be avoided by using large dynamic range sensors. Usually, large measuring range (with the same sensor resolution) sensors tend to be expensive. Instead of using an expensive large range sensor, we can use a low cost one with smaller measuring limit along with applying a saturation compensation algorithm.
Similarly, we can define S g,y , S g,z , y a,sat (saturation value of an accelerometer), S a,x , S a,y and S a,z . Denote δ g,x,k , δ g,y,k , δ g,z,k , δ a,x,k , δ a,y,k , δ a,z,k , the compensation values of gyroscopes and accelerations in ,, x y z axes at the time k , respectively. The compensated sensor value ( , gk y and , ak y ) is the sum of sensor output value and compensation value. For example, the compensated x axis gyroscope value is Let  be set of ,,  smoother is can be evaluated using the computed J value in Equation (15). This process (steps B and C) is repeated by changing  values. The algorithm finishes if the minimum value of J is found. We note that  optimization can be formulated as a constrained nonlinear optimization problem. Once the minimization process is done, we can compute the saturation compensated smoother values (step E). In this section, f is minimized with respect to .
 Two different methods are proposed. The first method directly minimize  with respect to all possible combination of .
 The second method uses the geometric structure of the sensor saturation.

Method 1: Direct Estimation of 
In the first method, the following optimization problem is solved: subject to: In Equation (20) Compute the saturation compensated smoother value using (4) and (7) quadratic optimization constrained nonlinear optimization A B C D E gyroscope data, the maximum angular velocity of a human knee varies from 213 to 1,087 s  [22]; in gait acceleration, the maximum foot acceleration is around 11.82 2 ms for walking case [23]. Therefore, the gyroscope data upper bound can be chosen as δ b = [1,087y g,sat 11.82y a,sat ] T . From Equations (16) and (17), it is easily to be seen that when ,0 ga yy , we have 0   because gyroscope compensated values ( , ga yy ) are larger than saturation value ( ,, , g sat a sat yy ). The initial value of  can be chosen as a set of 0 and varies in a range which satisfies the condition in Equation (20). With each set of  , one value of x is obtained. The minimum value of  which makes x minimizes the quadratic problem (14) is chosen.

Method 2: Estimation of  Using Geometric Form
In this method, the saturated sensor data region is approximated by a triangle (see Figure 4a) and quadratic function (see Figure 4b). Using three data before the saturation interval ( 2 The other compensation value is defined by: Using this method, sets of sensor's compensation values  can be generated by changing the intersection point's value from saturation value to maximum compensation value. In this case,  only depends on intersection point. Denote

 
, II I x t the intersection point. Once I is created, the rest compensation data will be obtained based on 1 l and 2 l lines using Equation (22). The optimization problem (19) will subject to  value which relates to the intersection point. With the same idea, a quadratic approximation can be used to generate the lost information (see Figure 3b. Firstly the maximum compensation value can be generated using the same procedure in triangle approximation process. After this step, a quadratic function () t

Sensor Saturation Compensation Algorithm with Multiple Zero Velocity Intervals
Section 4 introduced two methods to compensate sensor saturation in a standard movement situation which contains two zero velocity intervals at the beginning and ending of a moving interval. In this section, we apply the Section 4 methods to multiple zero velocity interval movements. A movement with multiple zero velocity intervals is displayed in Figure 5. In this movement, moving intervals are interposed by zero velocity intervals. An example of this movement can be found in gait analysis in [4,5]. In the characteristic human gait, one foot is assumed to be not moving when the the center of mass is put on the corresponding leg to move the other. Therefore, during the walking process, there are zero velocity intervals separated by moving intervals for each foot.
We divide the movement into segments based on zero velocity intervals so that between two zero velocity intervals there is data from one moving period (see Figure 5). If there are saturations in the movement, we can apply the compensation methods in Section 4 to each segment. The information of the last sampling in the prior segment will be used as the initial value for the following one. The initial error covariance for each segment can be estimated using a Kalman filter with a zero velocity measurement update for previous segment. Note that for the first segment the position and velocity errors are assumed to be zero due to the fact that the initial body coordinate frame coincides with the navigation coordinate frame. The initial covariance of the first segment is also assumed to be small. Repeating the smoothening processes until the last segment, the whole smoothed data is obtained.

Simulation and Experiments
In this section, one simulation result and two experimental results are given to verify the proposed algorithm. First, a simulation is done to verify the proposed method. The sensor is assumed to be located at the end of a bar, which is rotated along the body y axis. There is a sensor's saturation in the gyroscope y axis data as in Figure 6. The saturated data (green "-" line) is obtained with ,, 7.5 g y sat y  . This saturated data is used as an input data for the compensation algorithms in Section 4.  In this simulation,  only contains y axis compensation (δ g,y,k ). The smoothing results (forward-backward smoother [20] and the proposed methods) are given in Figure 7 and Table 1. In case of method 1, δ b is chosen as 5, so that the performance index J is minimized subject to 05  . Similarly, method 2 (both triangle and quadratic approximation) is applied. In Table 1, we can see that both method 1 and 2 produce significant improvements over the conventional smoothing (forward-backward smoother) result. Among the proposed methods, method 1 gives the best result and the second best is method 2 (quadratic approximation). We note that the method 1 requires more computation than method 2 since it needs more optimization variables. As for method 2, it is not conclusive whether quadratic approximation is better than triangle approximation since triangle approximation sometimes gives better results as seen in our later experiment (see Table 2). In Figure 8, the gyroscope sensor saturation compensation results by the proposed methods are given. In all three cases, it can be seen that the sensor saturation is well estimated.  Figure 7, the gyroscope sensor saturation compensation results by the proposed methods are given. In all three cases, it can be seen that the sensor saturation is well estimated.
In order to verify the robustness of the proposed algorithms to the saturation, a simulation is done by checking the position norm errors while ,, g y sat y value is changed. For example, if we set y g,y,sat =12 in Figure 6, there is no saturation in the sensor. On other hand, if we decrease y g,y,sat value, the sensor saturation increases. With the method 1, when the b  is fixed and the saturation value is changed, the accuracy of the proposed smoother may be affected (if the b  is chosen not large enough). As can be seen from Figure 9 where the b  is chosen as 2.5, the accuracy will decrease when the compensation value that is used to compensate the saturated data is larger than b  . However, this effect can be avoided since the saturation value is known for each sensor. Therefore, we can choose b  as a large number compared with the sensor's saturation value.  Figure 10 illustrates that the saturation is well compensated by method 2, even when the saturation value is changing. In general, compared with a forward-backward smoother, our proposed compensation smoother has better performance.
As mentioned above, how to choose b  may affect the accuracy. In some applications, the b  is known from doing statistical experiments. In the case where b  is unknown, we can choose an arbitrarily large value. This does not affect the result as shown in Figure 11. In Figure 11, when b  is chosen smaller than the compensation value that is used to compensate the saturated data, the error could be large. If we increase b  , the error decreases since a larger saturation can be compensated. However, it could lead to more computations for solving the optimization Equation (19).  In the first experiment, an object, which is attached with an IMU on the top, moved in a straight line of 0.95 m so that the x axis of the IMU coincides with direction of movement. In this case, there exists saturation in the sensor's accelerometer in the x axis ( , ax y ). The trajectory of the object is estimated by the forward-backward smoother and our proposed compensation smoother using two methods. Figure 12 shows that the saturated ,, a x k y data was compensated using the proposed compensation smoothers. The estimated position errors are given in Figure 13 (method 1 only, since the results of method 2 are similar) and Table 2. As can be seen from Figure 13, the proposed compensation smoother using method 1 gives a closer result to the true position (0.9121 m, equivalent to 4% error) than the forwardbackward smoother (0.8404 m, equivalent to 11.54% error). The proposed compensation smoother using method 2 also gives a good result of 0.9084 m (4.39% error) and 0.8976 m (5.52% error) for triangle and quadratic approximations, respectively.  Another experiment has been done to verify the compensation feasibility of the proposed algorithm. In this experiment, an IMU is attached at the tip of a digitizer (see Figure 14). The tip was moved in a curved line. The position data of the tip was recorded by the digitizer while the trajectory of the IMU is estimated by the proposed compensation and the forward-backward smoothers, respectively. The movement was made so that there is saturation in the ,, g z k y data. The estimated trajectories are compared with the true trajectory of the digitizer in Figure 15 (in millimeters). The distance between the start and stop point is used as an evaluation criterion. Based on this criterion, a table of distance errors of the estimated positions and the digitizer's data is formed, as shown in Table 3. Table 3 shows that method 1 gives a best accuracy (7.9 mm error) compared with other estimations. The forward-backward smoother has a worst estimation with a 50.4 mm error result. Moreover, the two approximation methods in method 2 provide similar results (11.8 and 10.3 mm errors).   In the last experiment, we verify the compensation smoother application in multiple zero velocity intervals movement (in Section 5). In this experiment, an IMU is attached on a human foot. The volunteer was asked to walk along a straight corridor. A pen is also attached on the volunteer's shoe to mark the steps' positions on the floor. The obtained data from IMU is used to estimate the trajectory of the foot. A comparison of proposed and forward-backward smoothers trajectories is given in Figure 16. The result shows that the last position error of proposed smoother trajectory is 0.9042 m, while it is 2.3128 m for the forward-backward smoother trajectory.

Conclusions
This paper has proposed some approaches to compensate for sensor saturation. Saturation is a common problem with IMU sensors in tracking a moving object. The lost data in the saturated parts could be important due to the accumulated errors. In the paper, the authors used a standard smoothing algorithm with zero velocity intervals to compensate sensor saturation. The considered motion includes a moving interval between two zero velocity intervals. Two methods were proposed. The first method directly estimates the saturation compensation while the second one uses a geometric form to estimate the saturation. The proposed smoothing algorithm can be applied in some motions which contain many moving intervals separated by zero velocity intervals. In this case, the motion is divided into segments based on zero-velocity intervals so that between two zero velocity intervals there is one moving interval. The saturation estimation algorithm is applied in each segment from the first to the last one. To verify the feasibility of the two methods, some experiments have been done. The experiments showed that the proposed smoothing algorithm can compensate the sensor saturation and provides a smaller error than a conventional smoother (forward-backward filter). In practical applications, the sensor saturation compensation methods proposed in this paper can be used to improve the accuracy of small dynamic range sensors instead of using a large dynamic range one which usually tends to be more expensive.