Inertial sensor-based smoother for gait analysis.

An off-line smoother algorithm is proposed to estimate foot motion using an inertial sensor unit (three-axis gyroscopes and accelerometers) attached to a shoe. The smoother gives more accurate foot motion estimation than filter-based algorithms by using all of the sensor data instead of using the current sensor data. The algorithm consists of two parts. In the first part, a Kalman filter is used to obtain initial foot motion estimation. In the second part, the error in the initial estimation is compensated using a smoother, where the problem is formulated in the quadratic optimization problem. An efficient solution of the quadratic optimization problem is given using the sparse structure. Through experiments, it is shown that the proposed algorithm can estimate foot motion more accurately than a filter-based algorithm with reasonable computation time. In particular, there is significant improvement in the foot motion estimation when the foot is moving off the floor: the z-axis position error squared sum (total time: 3.47 s) when the foot is in the air is 0.0807 m2 (Kalman filter) and 0.0020 m2 (the proposed smoother).

Introduction

Gait analysis is the study of human walking [1-3]. It is helpful in the medical diagnosis and treatment of those diseases that affect the ability to walk.

There are many devices that can be used to perform gait analysis [4]. Mechanical goniometers are used to measure hip, knee and ankle angles. Force plates on the floor or on a shoe (such as the wearable force plate from Tec Gihan Co., Kyoto, Japan) are used to measure the vertical force beneath the foot [5]. Optical motion trackers (such as Vicon) are used to measure the movement of a leg and a foot. These sensor systems can be used independently or together to give necessary gait parameters, such as stride length, walking speed and ground reaction forces.

Probably, an optical motion tracker is the most popular gait analysis tool [6], since it provides a very accurate description of foot and leg motion. The disadvantages of an optical motion tracker are its high cost and limited walking range. You can only increased walking range by using additional cameras, which can be quite expensive for a long walking range.

An alternative technology for the optical motion tracker is the use of inertial sensors [7,8]. Although the accuracy of inertial sensor-based estimation is generally worse than that of the optical motion tracker, the inertial sensor approach does not have a walking range limitation, and the cost is generally cheaper. If inertial sensors (gyroscopes and accelerometers) are attached to a leg, the angles of leg joints can be estimated (see [9-11] and commercial products from INSENCO and Xsens). These angle-based systems are extended to estimate whole human body motion [12-14]. Theoretically, the position of a leg also can be estimated using the inertial navigation algorithm [15]. However, the position estimation error quickly diverges when low cost inertial sensors are used. If inertial sensors are attached to a foot, the position (in addition to angle estimation) can be estimated using an inertial navigation algorithm. Unlike in the leg case, we know that the foot touches the ground almost periodically during walking. When the foot touches the ground, the velocity of the foot becomes temporarily zero, and during this zero velocity interval, the velocity errors of the inertial navigation algorithm can be reset. This zero velocity updating significantly reduces the position error [16].

Estimation of foot movement by attaching inertial sensors to the foot has been extensively studied in the context of a pedestrian navigation system [17,18]. In the filter-based foot motion estimation, the position error increases as soon as the foot is off the ground, since there are no measurements when the foot is in the air. The position error decreases once the foot touches the ground from the zero velocity updating. Thus, in the filter-based foot motion estimation, the motion estimation is not accurate when the foot is moving, but the final position (after the zero velocity updating) becomes accurate. In the pedestrian navigation system, the final position of the foot (that is, the final position of a person) is the main concern and, thus, the filter-based estimation can be used.

On the other hand, foot motion trajectory is the main concern in gait analysis. The filter estimated position can be improved using the smoother. The filtering algorithm uses the measurement data up to the current time to compute the current estimation. On the other hand, the smoother uses the whole measurement data (including the future data) to compute the current estimation. In [19], a smoother algorithm is proposed based on the forward-backward filter, where a smoothing algorithm is applied for each walking step. The position after the zero velocity updating is used as the initial value of a backward filter. The smoothed result for each step is combined rather heuristically to produce total walking estimation.

In this paper, the smoothing problem is formulated as a global optimization problem, where the whole walking motion is estimated in a single optimization problem. The smoother algorithm is formulated as a quadratic optimization problem (motivated by [20]), where its problem size depends on the total walking time. Since the problem size becomes very large even for a few minutes of walking data, the efficient solution is given utilizing the sparse structure of the problem.

The remainder of the paper is organized as follows. In Section 2, a basic material is introduced defining the coordinate systems and the main variables. In Section 3, a Kalman filter to estimate the foot motion is given. In Section 4, a smoother algorithm to estimate the foot motion is formulated in the quadratic optimization problem. In Section 5, a solution to the quadratic optimization problem is given using the sparse structure of the problem. In Section 6, experiment results are given to verify the proposed algorithm. Concluding remarks are given in Section 7.

Problem Formulation

Two-coordinate frames (the body coordinate frame and the navigation coordinate frame) are defined as in Figure 1. The three axes of the body coordinate frame coincide with the three axes of the inertial sensor unit. The z-axis of the navigation coordinate frame coincides with the local gravitational direction. The choice of the x-axis of the navigation coordinate frame can be chosen arbitrarily, since the walking direction is not important in the gait analysis. The notation [p] ([p]b) is used to denote that a vector p is represented in the navigation (body) coordinate frame.

Figure 1.

Overview of the foot inertial sensor unit.

The position of the foot is defined by [r] ∈ R3, which is the origin of the body coordinate frame expressed in the navigation coordinate frame. Similarly, the velocity of the foot is denoted by [v] ∈ R3. The attitude of the foot is represented using a quaternion q ∈ R4, which represents the rotation relationship between the navigation coordinate frame and the body coordinate frame. The directional cosine matrix corresponding to quaternion q is denoted by C(q) ∈ SO(3).

The basic dynamic equations [21,22] of q, r and v are given by: where ω = [ω ]′ is the angular velocity of the body coordinate frame with respect to the navigation coordinate frame. The symbol Ω(ω) is defined by:

The external acceleration (acceleration related to the movement, excluding the gravitational acceleration) expressed in the body coordinate frame is denoted by a ∈ R3.

The goal of this paper is to estimate q, r and v. The estimated values are denoted by q̂, r̂ and v̂, respectively The inertial navigation algorithm [19,23] is used to estimate q, r and v, where the equations in Equation (1) are integrated to compute q̂, r̂ and v̂. Accelerometer output y ∈ R3 and gyroscope output y ∈ R3 are related to ω and a as follows: where n ∈ R3 and n ∈ R3 are sensor noises. The vector g̃ is the local gravitational acceleration, where g̃ = [0 0 9.8] is used. The sensor biases are denoted by b ∈ R3 and b ∈ R3. It is assumed that the sensor biases are constant and unknown.

Let T be the sampling period of the sensor output. The sampled outputs of y and y are denoted by y = y(kT) and y = y(kT). The same notation is used for q, r and v. Let N be the final discrete time index of the sensor output, where the starting index is one. The discrete sensor noise n and n are assumed to be zero mean white Gaussian noises, whose covariances are given by:

Due to high sensor noise level of low cost inertial sensors and unknown biases, an inertial navigation algorithm without external reference measurement can estimate foot motion only for a few seconds. To overcome this problem, a zero velocity updating method is used in an inertial sensor-based foot motion estimation. During walking, the foot touches the ground almost periodically. When the foot is on the ground, we know that the velocity of the foot is zero. Thus, whenever the foot touches the ground, the velocity error in the inertial navigation algorithm can be reset, which significantly reduces foot motion estimation errors. To use this zero velocity updating, zero velocity intervals (that is, when the foot is on the ground) must be detected. Many zero velocity detection algorithms have been proposed [16,24].

In this paper, we use a simple zero velocity detection algorithm. Let Z be a set of all discrete time indices belonging to the zero velocity intervals. The discrete time k belongs to Z if the following conditions are satisfied: where N and N are even number integers.

Kalman Filter

In this section, q, r, v, b and b are estimated using a Kalman filter. The estimated values are denoted by q̂, r̂, v̂, b̂ and b̂. As in [19,23], q ∼ b are not directly estimated in the Kalman filter. They are first estimated from the integration equation of (1), and then, the errors in the integrated values are estimated in the Kalman filter. These estimated values will be used in Section 4 to derive a smoothed estimation.

Initialization

Since only the relative position is important in the gait analysis, the initial position is assumed to be zero: that is, r̂,1 = 03×1, where “1” denotes the initial discrete time index. Assuming that a person is not moving at the initial time, the initial attitude is computed using the triad algorithm [25,26]. The following two vector relations are used in the triad algorithm:

The pitch and roll angles are determined from the first vector relation, and the yaw angle is determined from the second vector relation. Since the heading of walking is not important in the gait analysis, the yaw angle can be chosen arbitrarily.

The initial velocity is chosen to be zero: that is, v̂,1 = 03×1. If there is no particular information on b and b, the initial estimates are chosen to be zero (b̂ = b̂ = 03×1).

Integration of Estimated Values

Suppose we have computed q̂, r̂, ⋯ , b̂ at the discrete time k. The estimated values at the discrete time k + 1 are firstly computed using Equation (1), where ω and a are replaced by y − b̂ and y − C(q̂)g̃ − b̂:

For later use, we define a function f, which represents a numerical integration of Equation (6) from kT to (k + 1)T:

We assume that the variables not used as arguments of f are embedded in f: for example, the ω, C(q) and g̃ terms are assumed to be embedded in the function f.

Kalman Filter

Due to noise terms (n and n) and bias estimation error (b − b̂ and b − b̂), there are errors in the numerically-integrated values of Equation (7). These errors are estimated using a Kalman filter. Let X be defined by: where ⊗ is the quaternion multiplication and q* denotes the quaternion conjugate of a quaternion q [21]. The elements of X represent the estimation errors in q̂ ∼ b̂. In defining X, it is assumed that the quaternion errors in q̂ are small. Thus, (q̂)* ⊗ q is approximated by:

With this assumption, the error in q̂ ∈ R4 is represented by three-dimensional vector q̄ instead of a four-dimensional vector. See [27] for the details about the multiplicative error notation of a quaternion. The initial estimate is set to X̂,1 = 015×1, and the initial error covariance P1 is given by:

The initial attitude error covariance P is given by the algorithm in [25], where R and 10−11I3 are used as measurement noise covariances of the first and second equations in Equation (5). The small covariance value for the second equation is to make P nonsingular, since

is used in the smoothing algorithm. The remaining covariances P ∼ P can be considered as design parameters, whose values are given in Table 1.

Table 1.

Parameters of the proposed smoother and Kalman filter.

Parameter	Symbol	Values	Related Equation
sampling period	T	0.01
zero velocity detection parameters	Bg, Ng, Ba, Na	0.2, 10, 0.4, 10	(4)
sensor noise error covariance	Rg, Ra	0.001I3, 0.01I3
measurement noise covariance	Rv, Rr	0.0001I3,0.0001	(15)
filter and smoother parameter	Pr,init, Pv,init	0.0001I3, 0.0001I3	(9), (27)
filter and smoother parameter	Pbg,init, Pba,init	0.000001I3, 0.000001I3	(9), (27)

The dynamic equation of X is given by: where:

For a vector p = [p1 p2 p3]′ ∈ R3, [p×] is defined by:

Note that A is a lower triangular matrix, and thus, e is also a lower triangular matrix having the structure in Equation (11) [28].

The process noise w is the zero mean white Gaussian noise with: where the 10−6 value in Q is added to make Q nonsingular, since is used in the smoother. Derivation of A and Q can be obtained by slightly modifying the results in [19,23]. The computational load of e can be reduced by using the following approximation:

Remark 1

How the estimation error in k-th discrete time evolves after the integration Equation (7) is represented in Equation (10). For example, (10) states that the error in r̂+1 comes from two sources. One is from the estimation error in k-th time (before the integration), and the other is from the process of numerical integration. The latter is not from the numerical integration error (it is ignored), but from the noise terms and bias estimation errors.

In an inertial sensor-only foot motion estimation, there is no external measurement. The fact that the velocity is zero during the zero velocity intervals is used as a fictitious measurement.

During the zero velocity interval, the following equation can be used as a measurement equation: where n is a fictitious measurement noise representing the confidence of zero velocity interval detection in Equation (4). The z-axis value of r̂ is almost the same when the foot is on the ground, assuming that a person is walking on a flat floor. Thus, the following measurement equation can also be used during the zero velocity interval.

where n is a fictitious measurement noise.

We assume that n and n are zero mean white Gaussian noises, where covariances are given by:

A standard discrete time Kalman filter is applied to Equations (10), (13) and (14). If k ∈ Z, we can compute X̂ from the Kalman filter, and these values are used to update q̂ ∼ b∈.

The Kalman filter algorithm of position estimation (r̂) is summarized in the following for k1 ≤ k ≤ k2 assuming that k1 ∈ Z, k2 ∈ Z and k ∉ Z for k1 < k < k2. The algorithm for other values (q̂, v̂ and b̂) are similar.

discrete time k = k1 ∈ Z

r̂ is given from the previous step

r̄̂ = 0

discrete time k1 < k < k2

compute r̂ using Equation (7)

r̄̂−1 = 0 from Kalman filter time update Equation (10) with r̄̂−1 = 0

discrete time k = k2 ∈ Z

compute r̂ using Equation (7)

compute r̄̂ from the Kalman filter measurement update Equations (13) and (14)

r̂ ⟸ r̂ + r̄̂

r̄̂−1 = 0

The result in this section is not a new result, and similar results can be found in [19].

Smoothing Algorithm

In this section, the Kalman filter-estimated values (q̂ ∼ b̂) in Section 3 are improved using the smoother, where the errors in q̂, r̂, v̂ and b̂ are compensated for. To do that, the estimation error is defined as follows. where b̂ and b̂ are b̂ and b̂ (the final values) in the Kalman filter, which can be considered as the most accurate estimation of b and b. Since b and b are constant, b̄ and b̄ are not time dependent, while q̄, r̄ and v̄ are time dependent. We note that the estimation errors in Equation (8) represent the errors before the measurement update, while the errors in Equation (16) represent the errors after the measurement update.

Let X be defined by:

Equations for time-dependent variables (q̄, r̄ and v̄) are given by: where Φ and Ψ are defined in Equation (11) and:

Derivation of Equation (18) is now explained. From Equations (16), (19) and (20), we have:

The second and the third elements of Equation (21) are straightforward. The first element of Equation (21) is derived from the following: where the following approximation is used:

From Equation (21), ξ + X+1 denotes the estimation errors of f in Equation (20). Thus, Equation (18) represents how the estimation error evolves after the integration Equation (20). Thus, Equation (18) is essentially the same as Equation (10). In the exact same methods in the case of Equation (10), Equation (21) can be obtained by slightly modifying the results in [19,23].

The smoothing problem is formulated as a quadratic optimization problem using the method in [20]. Let an optimization variable X̃ be defined by:

The smoothing problem can be formulated as the following quadratic optimization problem: where:

Initial information on b and b are reflected in b ∈ R3 and b ∈ R3. For example, b can be computed if a sensor unit is put on the ground without moving for 1∼2 min. If there is no information on b and b, we set b and b to zero.

Inserting Equation (18) into Equation (27), we have:

It is not difficult to see that Equation (28) is a quadratic function of X̃ ∈ R9×+6: where M1, M2 and M3 can be computed from Equation (28). M1 and M2 are given in the Appendix (M3 is omitted, since it is irrelevant in the optimization). The minimizing solution of Equation (29) can be computed by solving the following linear equation: where X̃* is the minimizing solution.

Once the minimizing solution X̃* is computed, the smoother estimated values are computed using Equation (16). For example, the position estimate is computed as follows: where r̂ is the smoother estimate of r and is the 9(k − 1) + 1 ∼ 9k rows of X̃*.

The result in this section is a new result. The main contribution of this section is that the foot motion Kalman filter error is compensated for using the quadratic optimization problem.

Smoothing Optimization Solution

Although the smoothing problem boils down to a simple linear Equation (30)M1 ∈ R(9 becomes very large when N is large. Suppose we solve Equation (30) using the Cholesky decomposition. Let L be the Cholesky triangle (lower triangular matrix) satisfying:

Once L is computed, X̃* can be computed by the following equation:

The computational load of Cholesky decomposition Equation (31) is (9N + 6)3/6 flops, and the load for each linear equation in Equation (32) is (9N + 6)2/6 flops [29]. Suppose T = 0.01 and N = 12,000 (2-min walk data), then the computational load of the Cholesky decomposition alone is approximately 2 × 1014 flops, which is practically difficult to compute (due to the very long computation time and large memory requirement).

Fortunately, the matrix M1 has a sparse structure, and the computational load can be reduced using the sparse structure.

We partition M1 and M2 as follows:

We also partition X̃ in Equation (29) accordingly:

From the definition of M11 in the Appendix, we can see that M11 is a banded matrix [29] with the bandwidth p = 17: where * denotes the arbitrary 9 × 9 block matrix. All unspecified entries are zero 9 × 9 block matrices.

Equation (30) can be rewritten using the partitions in Equations (33) and (34):

Equation (36) can be transformed into the following equation using the Schur complement [30]:

The Equation (37) can be solved in the following sequence [30]. First, find the Cholesky decomposition of banded matrix M11: where L11 ∈ R99 is a lower triangular matrix. Find U1 ∈ R9×1 and U2 ∈ R9×1 from the following linear equation:

Inserting Equation (39) into Equation (37), we have:

After computing X̃2 in Equation (40), we compute X̃1 from the following equations:

When solving Equations (38)–(41), the algorithms for banded matrices in [29] can be used for efficient computation and memory usage. The computational load of the Cholesky decomposition of the banded matrix Equation (38) is (9N(p2 + 3p))/2 flops (p = 17 is the bandwidth of M11, and 9N ≫ p is assumed) and 9N square root computation. The computational load of each equation in Equations (39) and (41) is 9N(p + 1) flops. Equation (40) is a small dimensional problem (X̃2 ∈ R6).

Compared with the standard method, the solution using the banded structure of M11 requires much smaller computation and memory. For example, the Cholesky decomposition requires approximately 3.6 × 107 flops and 9N square root computation in the case of N = 12, 000. Compared with 2 × 1014 flops in the standard computation case, its computational requirement is significantly small.

The main contribution of this section is to derive an efficient solution to the quadratic optimization problem and to investigate the quantitative computational requirement.

Experiment

Two experiments are performed to verify the proposed algorithm. In the experiments, an inertial sensor unit (Xsens MTi sensor unit) is attached on the foot, as shown in Figure 2. The parameters used in both experiments are given in Table 1.

Figure 2.

Experimental setup for the first experiment.

In the first experiment, a person walked three steps, while the movement of the foot is tracked using an optical motion tracking system (Optitrack six Flex 13 camera system). The experimental setup is given in Figure 2. Six optical markers are used in the optical tracking system, and the center position of six markers is used as the true foot position. Since the coordinate systems of the optical tracker and the inertial sensor unit are different, the optical tracker output is rotated and translated.

The estimated position (Kalman filter and the proposed smoother) with the flat floor assumption (see Equation (14)) is given in Figure 3: the left figure is the z-axis movement, and the right figure is the xy-axis position. In the left figure, the zero velocity intervals are represented by the “*” symbol.

Figure 3.

Estimated position with the flat floor assumption: z-axis position (Left) and xy-axis position (Right).

Assuming that the motion tracker value is the true position value, we can see that the Kalman filter gives a good final position estimation (around 5 s). However, the position estimated during walking is not accurate; for example, see the values between 2 and 3 s. Thus, the Kalman filter result may be good enough for the pedestrian navigation system, where the final position is important. However, the Kalman result is not good enough for gait analysis, where it is important to know how the foot has moved in addition to the final position. The final estimated position of the smoother is almost similar to that of the Kalman filter. However, we can see that there is significant improvement in the middle foot position estimation.

The z-axis estimation error is compared quantitatively in Table 2. The z-axis estimation error e is defined by: where r̂ is either r̂ = r̂ or r̂ = r̂.

Table 2.

Estimated z position error squared sum with the flat floor assumption.

Algorithms	E_1=∑_k∈Z_mez,k2	E_2=∑_k∉Z_mez,k2	E1 + E2
	(error squared sum when the foot is on the ground)	(error squared sum when the foot is in the air)
Kalman filter	0.0008	0.0807	0.0815
Proposed smoother	0.0009	0.0020	0.0029

The z-axis estimation error squared sum is given in Table 2. The estimation error of both algorithms during the zero velocity interval (k ∈ Z) are almost the same, since the estimated values becomes zero with the flat floor assumption. The estimation error when the foot is in the air are different; the proposed smoother gives significantly smaller estimation error: 0.0020 in the case of the proposed smoother and 0.0807 in the case of the Kalman filter. Furthermore, the maximum z-axis position errors for the Kalman filter and the proposed smoother are 4 cm and 0.5 cm, respectively. This improved estimation is due to the fact the smoother uses all of the available information.

The position estimation by the Kalman filter and the proposed smoother without the flat floor assumption are given in Figure 4. Without the flat floor assumption, we can see that the z-axis error increases since the z-axis position error is not corrected using Equation (14).

Figure 4.

Estimated position without the flat floor assumption: z-axis position (Left) and xy-axis position (Right).

Due to the viewpoint constraint of the optical tracker, the walking distance is rather limited in the first experiment. In the second experiment, a person walked 50 m along a straight line on the flat floor. In this experiment, there is no ground truth value, except the z-axis value, which should be almost constant during the zero velocity intervals, since a person walked on a flat floor.

The Kalman filter estimated and smoother estimated positions are given in Figure 5, both with the flat floor assumption. In the Kalman filter result, the estimated xy trajectory is slightly curved. The smoother estimated xy trajectory is a little bit more straight. The true trajectory should be almost a straight line, since the person walked along a straight line.

Figure 5.

Estimated position (50-m walk): z-axis position (Left) and xy-axis position (Right).

Similar to the first experiment, there is no significant improvement in the smoother position estimation for the position where the foot is on the floor (that is, after the zero velocity updating). The smoother, however, gives more accurate foot motion estimation while the foot is moving in the air, which is important in the gait analysis. This can be seen in the z-axis graphs: the Kalman filter-estimated z-axis trajectory is not accurate, since the z-axis value cannot be negative (which means the foot is beneath the floor).

The 50-m walk is repeated by four subjects: five times for each person, thus a total of 20 experiments. The total walk length estimation is given in Table 3. The average estimated walk length is 48.99 m with the standard deviation being 0.26 m. From the small standard deviation, the estimation result does not seem to depend on the subjects.

Table 3.

Four persons' 50-m walk experiment results.

Subject Number	Final xy Position Estimated Value (m)	Average (m)
1	48.97, 48.97, 49.16, 49.06, 49.02	49.04
2	49.36, 49.29, 49.21, 49.37, 49.14	49.27
3	48.52, 48.98, 48.46, 48.69, 48.45	48.62
4	48.98, 49.00, 48.92, 49.08, 49.06	49.01

We note that the xy-axis position error increases as the walk length increases. This is unavoidable, since the motion is estimated using the inertial navigation algorithm without any external reference. The z-axis position error increment is limited with the flat floor assumption.

The computation time of the proposed smoother was 8.7 s (MATLAB running on Intel Core i5-2310 CPU Windows PC, 2.9 GHz), where the walk time is 58.92 s (N = 5892).

Conclusions

In this paper, a smoother algorithm is proposed for gait analysis to obtain more accurate foot motion estimation than the filter-based estimation. In the smoother algorithm, the motion estimation problem is formulated as a quadratic optimization problem, and an efficient computational algorithm is given.

In the first experiment, we have compared the results of the Kalman filter and the proposed smoother. The Kalman filter-based estimation is quite accurate after the zero velocity update. However, foot motion estimation is not accurate when the foot is in the air. We have shown that the motion estimation is significantly improved for the foot motion in the air if the smoother is used: the z-axis position error squared sum when the foot is in the air is 0.0807 (Kalman filter) and 0.0020 (the proposed smoother). Thus, the proposed smoother is a useful algorithm for gait analysis, where the foot motion estimation is important.

In the second experiment, we investigated a long walk foot motion estimation (50 m). There is not much difference in the final position estimation of the Kalman filter and the proposed smoother. The 50-m walk was tested using four subjects, and we found that the estimation does not depend on the subject. This is not surprising, since there are no person-dependent parameters in the algorithm.

Since the proposed smoother algorithm is more complicated than the Kalman filter, the computational time is slow. In the second experiment, the computation time of the proposed smoother is about 8.7 s to process about one minute of walk data (the sampling period is 0.01 s), while it took 0.8 s for the Kalman filter. Thus, the estimation performance is improved with the additional computation time.