Use of ant colony optimization and the Kalman filter to deduce thigh dip angle via acceleration and angular velocity sensing

Abstract

Background:

Acceleration and angular velocity sensors are commonly used in the measurement of gait parameters. A representative application calculates the limb segment dip angle. The rotation angle is typically deduced by a conventional Kalman filter, which includes the use of two empirically derived parameters. We improved this conventional method by introducing colony algorithm to find the optimal parameter combination instead of empirically assignment.

Method:

To achieve optimal results, a servo motor with an inertial measurement unit was used to simulate human limb segment motion according to programmed rotation angles that was employed as the ground truth. To minimize the bias between the ground truth and the calculated result, the ant colony algorithm was employed to obtain the optimal Kalman filter parameter combination in two-dimensional space.

Results:

By the motor experiment, the sum of the angle squared error was only 1.9305 rad², much better than the 6.7723 rad² error by the conventional method. The optimal parameter combination obtained was then used in a human experiment involving a basketball player. The frames from video of a whole gait cycle period were all showed with a corresponding deduced thigh dip angle curve diagram.

Conclusion:

The colony algorithm for parameters optimization results in less angle errors deduced by Kalman filter using the data from inertial measurement unit. The subject experiment verified the feasibility and performance of this method.

Introduction

Since the 1970s, gait analysis has attracted considerable attention from researchers and clinicians. In human gait analysis, the measurement of lower extremity joint angles is very important and is useful in applications such as rehabilitation exercise, sports training, and entertainment. An inertial measurement unit (IMU) is the most commonly employed tool for estimating and analyzing kinematic movement. An IMU comprises accelerometers, gyroscopes, and magnetometers for measuring the acceleration, angular velocity, and the magnetic field vector in their local coordinate system, respectively. The orientation of the sensors can be estimated using appropriate algorithms, which can facilitate gait analysis based on an IMU mounted on a human body. An early method fused three accelerometers and a single gyroscope to monitor human posture and walking velocity.^1,2 To estimate the human posture in the sagittal plane, one of the three accelerometers was mounted on the subject’s trunk, thigh, and shank, respectively, thereby measuring the angle of the gravitational direction of each anatomical component. To estimate walking velocity, a gyroscope was attached to the thigh, and changes in the angle of the sagittal plane caused by walking were obtained by integrating the sensor signals. The walking velocity was then calculated according to changes in the thigh angle, the period of a single step cycle, and the subject’s leg length. Pappas et al. ³ proposed a reliable gait phase detection system to measure the angular velocity of the foot, and the forces exerted by the foot were assessed through a gyroscope and three force-sensitive resistors. Roetenberg et al. ⁴ used gyroscope, accelerometer, and magnetometer signals to estimate the orientations of various segments of a human body. Santhiranayagam et al. ⁵ mounted inertial sensors on a human foot to estimate gait parameters for evaluating their suitability in clinical applications. To obtain more information, such as for joint angles, additional sensors can be mounted on leg segments. Other researchers have employed a body-mounted gyroscope to calculate body segment orientations from the angular velocity data. ⁶ Another researcher employed three IMU sensors to track and locate human velocity. ⁷

In the 1980s, video camera systems were developed to track human movements, and, since then, this method has been an important method for gait analysis.^8,9 Gait movement can be analyzed using an optical tracking system, for example, the Vicon motion analysis system.^10,11 Andrews employed optoelectronic technology to measure the position of markers mounted to body segments.^12–14 However, this technology is expensive and is generally limited to indoor use, such as in laboratories or clinical environments. Moreover, all cameras employed in the approach require precise assembly. ¹⁵ However, the development of micro-electromechanical system (MEMS) technology has fostered a series of lightweight, low-power inertial sensors. ¹⁶ IMUs are based on this technology. Compared to video camera technology, MEMS technology is a better choice under particular conditions owing to its miniaturization, and the developed sensors are lightweight, portable, less expensive, safe, and much more easily mounted to the human body. ¹⁷ Wearable sensing systems have advantages over other systems, such as the capacity for conducting experiments outdoors to monitor daily movements. ¹⁸ One of the inconveniences is that the employment of a large number of sensors at one time may become overly burdensome for subjects. ¹⁹ Moreover, the sensors are subject to output errors mainly caused by noise, axis misalignment, and cross-axis sensitivity, particularly when the sensors are improperly handled.^20,21 The main output error of gyroscopes is associated with angle wander, which results from angle random walk and the progressive dynamic variation of the bias. An accelerometer produces noise in its output. Ferromagnetic materials, such as iron and other magnetic materials, near a magnetometer sensor will affect the direction and density of the local earth magnetic field and thus negatively affect the orientation estimates. ²²

Most existing studies have employed the Kalman filter to estimate the angle due to its better response.^23–26 However, the accuracy of a Kalman filter algorithm is highly dependent on the appropriate tuning of its parameters, and a balance suitable for both low- and high-intensity movements is difficult to achieve. In some studies, the parameters are empirically determined and achieve a tuning that may not be applicable to general conditions. Olivares and Górriz have proposed other filters such as least-mean squares (LMS) and recursive-least squares (RLS) filters. Miezal et al. ¹⁰ proposed an algorithm based on quaternions to calculate body segment orientation, which was converted from the orientation estimations of the sensor units through a rotation matrix obtained from a calibration trail. The general regression neural network (GRNN) algorithm has also been applied to estimate ankle, knee, and hip joint angles with good accuracy. However, to acquire this functionality, this method must first be trained for individual settings prior to measuring the parameters. Therefore, among all the approaches considered, we selected the Kalman filter to eliminate the drift. A common method employed to calculate the angle is to double integrate the detected angular acceleration, although this method results in significant drift.^27,28 To minimize this drift, a number of approaches have been proposed in the past decades. For example, Morris ²⁹ equated the signals obtained at the beginning and end of each walking cycle. Tong and Granat ³⁰ applied a low-cut, high-pass filter on the shank and thigh inclination angle signals. However, these methods cannot estimate the data correctly when the object is static or exhibits a low moving frequency. Thus, neither method can be used to process real-time signals. Another method for estimating joint angles from measured accelerations is to estimate shank segment orientation in the sagittal plane during walking and subsequently calculate the thigh inclination angle by adding the two values of the shank and knee angles. The results are acceptable only if the segment acceleration is small compared to that of gravity. ³¹

We utilized a six-axis IMU to obtain the original acceleration and angular velocity signals of the low limb segment. The Kalman filter was used to fuse the original data and calculate the segment angle. An experiment employing a standard subject was designed to verify the proposed algorithm. Simulated noise disturbance was added to the signal. Swarm intelligence algorithms have been widely used for parameter searching. ³² In this study, the ant colony algorithm was applied to obtain optimal values for two of the Kalman filter parameters according to ground truth criteria. The experimental results obtained for detecting the dip angle of an adult male thigh demonstrated a very high accuracy. This angle calculation algorithm is very fast and therefore suitable for real-time processing. The proposed method has good prospects in rehabilitation applications.

The key of this work is how to find the optimal parameters in the mathematical Kalman filter model. Most of the literatures have omitted the method of how to calculate the number of these parameters, while we have creatively utilized the ant colony algorithm to obtain the optimal parameters. A step motor was used as a standard subject to verify the angle deduction algorithm. A university basketball player was selected as the subject to show the result of thigh dip angle deduction method.

Hardware description

For this study of the human lower limb posture and motion during walking, we employed IMUs to measure the acceleration and angular velocity of thighs, shanks, and feet during walking. An MPU6050 six-axis IMU with a three-axis accelerometer and a three-axis gyroscope was employed in this study. The IIC protocol was employed for communications owing to its convenience in the data-acquisition experiment. Given that the data-processing system must accommodate multi-channel data acquisition, processing capacity, and real-time requirements, in addition to a secure digital (SD) card and a serial peripheral interface (SPI) communication protocol, we employed an ARM Cortex-M3 architecture STM32F103 microcontroller unit (MCU) as the main controller of the overall data-acquisition system to store data collected from IMUs. For the convenience of data acquisition, an appropriate data-storage method is important. The system stores the data provided by the IMU using an SD card storage method rather than a computer. The data are stored in the form of an Excel spreadsheet, in which each column represents a single axis of data. This type of stored file is easily read on a computer for subsequent offline data analysis.

Each data-acquisition experiment generated an Excel spreadsheet. The name of the spreadsheet comprises the date and the current number of the experiment on that day. Use of a real-time clock updates the date automatically, and the experiment number is automatically reset to zero when the date changes. The name of the current experiment is displayed on a flat-panel display for recording. The system employs a lithium polymer battery, which would be damaged if the voltage were lower than 2.85 V, and this necessitates the inclusion of a low-voltage alarm, such that the alarm is engaged when the battery voltage is less than 2.9 V as a reminder to replace the battery. The data-acquisition system detects the effectiveness of data acquisition each time after signal sampling. The system also checks the reliability of the hardware connection from the recorded file. If any IMUs are not connected, the system will sound an alarm and display an indication of the specific disconnected IMU(s). Experiments show that the performance of the data-acquisition system is reliable and stable. A schematic of the system is given in Figure 1.

Figure 1.

Hardware structure scheme.

Materials and methods

Angle deduction algorithm based on the Kalman filter

The quaternion coordinate system was used to represent an object’s attitude and changes in the three-dimensional (3D) angle. Unlike the Euler equations, gimbal lock does not occur when employing quaternions, which is convenient for solving the attitude. In this study, each limb segment was equipped with a single IMU with the goal of measuring the dip angle. After calculating the dip angles by updating the quaternion, in conjunction with the angular velocity data, the integration error was corrected by the Kalman filter using the angles calculated from the acceleration outputs.

Quaternion calculation

3D rotation using quaternions $\mathbf{q}$ can be expressed as follows

$$\mathbf{q}=\cos \frac{\mathrm{\phi }}{2}+\mathbf{n}\sin \frac{\mathrm{\phi }}{2}$$

(1)

$$\mathbf{r}\prime =\mathbf{qr}{\mathbf{q}}^{*}$$

(2)

$${\mathbf{q}}^{\mathbf{*}}=\cos \frac{\mathrm{\phi }}{2}-\mathbf{n}\sin \frac{\mathrm{\phi }}{2}$$

(3)

where $\mathrm{\phi }$ is the rotation of angle around unit vector $\mathbf{n}$ , ${\mathbf{q}}^{*}$ the conjugated quaternion of $\mathbf{q}$ , and $\mathbf{r}\mathbf{\prime }$ the vector after applying rotation to ordinary vector $\mathbf{r}$ .

The Runge–Kutta method was used to solve the following quaternion differential equation for updating the quaternion at each time t

$$\stackrel{·}{\mathbf{q}}={\mathbf{\Omega }}_{\mathbf{b}}·\mathbf{q}$$

(4)

$${\mathbf{K}}_1={\mathbf{\Omega }}_{\mathbf{b}}(t)·\mathbf{q}(t)$$

(5)

$$\mathbf{Y}=\mathbf{q}(t)+\Delta t·{\mathbf{K}}_1$$

(6)

$${\mathbf{K}}_2={\mathbf{\Omega }}_{\mathbf{b}}(t+\Delta t)\mathbf{Y}$$

(7)

$$\mathbf{q}(t+\Delta t)=\mathbf{q}(t)+\frac{\Delta t}{2}({\mathbf{K}}_1+{\mathbf{K}}_2)$$

(8)

$${\mathbf{\Omega }}_{\mathbf{b}}=\frac{1}{2}\left[\begin{array}{cccc}\hfill 0\hfill & \hfill -{\mathrm{\omega }}_x\hfill & \hfill -{\mathrm{\omega }}_y\hfill & \hfill -{\mathrm{\omega }}_z\hfill \\ \multicolumn{1}{c}{{\mathrm{\omega }}_x}& \hfill 0\hfill & \hfill {\mathrm{\omega }}_z\hfill & \hfill -{\mathrm{\omega }}_y\hfill \\ \multicolumn{1}{c}{{\mathrm{\omega }}_y}& \hfill -{\mathrm{\omega }}_z\hfill & \hfill 0\hfill & \hfill {\mathrm{\omega }}_x\hfill \\ \multicolumn{1}{c}{{\mathrm{\omega }}_z}& \hfill {\mathrm{\omega }}_y\hfill & \hfill -{\mathrm{\omega }}_x\hfill & \hfill 0\hfill \end{array}\right]$$

(9)

where $\Delta t$ is the sampling period, and the $\mathrm{\omega }$ terms are the components of the angular velocity vector $\mathbf{\omega }$ calculated from the gyroscope data. The value of $\mathrm{\phi }$ relative to the initial position can be obtained by equation (1), and the initial orientation of the sensor unit was estimated using acceleration data obtained from the human subject in an upright standing posture.

Sensor-to-segment mounting

As expressed before, the calibration posture was needed in this experiment. The experimenter was asked to stand upright for 5 s before walking. During that time, we needed to make sure the lower limb segment was straight and parallel to the direction of gravity. The relative position of the lower limb segment in the local sensor coordinate system, denoted by the vector ${\mathbf{v}}_s$ , can be represented by the acceleration data as

$${\mathbf{v}}_s=(\mathrm{ac}{c}_x,\mathrm{ac}{c}_y,\mathrm{ac}{c}_z)$$

(10)

The angle $\alpha$ between ${\mathbf{v}}_s$ and $\mathbf{n}$ can be calculated as

$$\alpha =<{\mathbf{v}}_s,\mathbf{n}>$$

(11)

Finally, the dip angle of the lower limb segment obtained from the angular velocity data, denoted ${\theta }_g$ , can be given as follows

$${\theta }_g=\mathrm{\phi }·\sin \alpha$$

(12)

Segment dip angle from acceleration data

The value of ${\mathbf{v}}_s$ is fixed, and the direction of gravity is obtained from the acceleration data of the initial segment orientation. The dip angle of the lower limb segment obtained from the acceleration data, denoted ${\theta }_a$ , can be given as follows

$${\theta }_a=<{\mathbf{v}}_s,{\mathbf{v}}_a>$$

(13)

$${\mathbf{v}}_a=(\mathrm{ac}{c}_x(t),\mathrm{ac}{c}_y(t),\mathrm{ac}{c}_z(t))$$

(14)

where ${\mathbf{v}}_a$ is the acceleration data obtained at time t. After obtaining both ${\theta }_g$ and ${\theta }_a$ , the Kalman filter was employed to fuse the two angles.

Application of Kalman filter

The Kalman filter uses a system’s dynamics model, known control inputs to that system, and multiple sequential measurements to form an estimate of the system’s varying quantities that is better than the estimate obtained using only one measurement alone. It is a common sensor-fusion and data-fusion algorithm. The Kalman filter deals effectively with the uncertainty due to noisy sensor data, and to some extent with random external factors as well.

The value of ${\theta }_g$ is very precise on short timescales, but, because even highly accurate gyroscopic sensors do not have a bias of exactly zero, the data are always subject to drift with time. The accelerometer data do not drift, so the value of ${\theta }_a$ would be expected to be more reliable and precise than ${\theta }_g$ over long timescales. However, ${\theta }_a$ is readily affected by accelerometer noise, which makes it less reliable under conditions of rapidly changing acceleration. To correct integration error, the Kalman filter is employed to fuse the two angles. Figure 2 presents a block diagram of the algorithm.

Figure 2.

Block diagram of angle calculation method using Kalman filter.

The state equation is represented by the error of the joint angle measured by a gyroscope, $\Delta \theta$ , and the bias offset of a gyroscope output, $\Delta b$ , as

$$\left[\begin{array}{c}\hfill \Delta {\theta }_{k+1}\hfill \\ \multicolumn{1}{c}{\Delta b_{k+1}}\end{array}\right]=\left[\begin{array}{cc}\hfill 1\hfill & \hfill \Delta t\hfill \\ \multicolumn{1}{c}{0}& \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill \Delta {\theta }_k\hfill \\ \multicolumn{1}{c}{\Delta b_k}\end{array}\right]+\left[\begin{array}{c}\hfill \Delta t\hfill \\ \multicolumn{1}{c}{1}\end{array}\right]w$$

(15)

where $w$ is the gyroscope measurement error. The observation signal is the difference between the gyroscope and accelerometer angles, $\Delta y$ , which is given as

$$\Delta y_k=\left[\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill \Delta {\theta }_k\hfill \\ \multicolumn{1}{c}{\Delta b_k}\end{array}\right]+v$$

(16)

where $v$ is the accelerometer measurement error. Based on this state-space model, the Kalman filter repeats corrections

$$\left[\begin{array}{c}\hfill \Delta {\hat{\theta }}_k\hfill \\ \multicolumn{1}{c}{\Delta {\hat{b}}_k}\end{array}\right]=\left[\begin{array}{c}\hfill \Delta {\hat{\theta }}_k^{-}\hfill \\ \multicolumn{1}{c}{\Delta {\hat{b}}_k^{-}}\end{array}\right]+\left[\begin{array}{c}\hfill K_1\hfill \\ \multicolumn{1}{c}{K_2}\end{array}\right](\Delta y_k-\Delta {\hat{\theta }}_k^{-})$$

(17)

and predictions

$$\left[\begin{array}{c}\hfill \Delta {\hat{\theta }}_{k+1}^{-}\hfill \\ \multicolumn{1}{c}{\Delta {\hat{b}}_{k+1}^{-}}\end{array}\right]=\left[\begin{array}{cc}\hfill 1\hfill & \hfill \Delta t\hfill \\ \multicolumn{1}{c}{0}& \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill \Delta {\hat{\theta }}_k\hfill \\ \multicolumn{1}{c}{\Delta {\hat{b}}_k}\end{array}\right]$$

(18)

where $K_1$ and $K_2$ are the Kalman gain for $\Delta \theta$ and $\Delta b$ , respectively, and the estimated and predicted values for $\Delta \theta$ and $\Delta b$ are represented by $\Delta \hat{\theta }$ and $\Delta {\hat{\theta }}_k^{-}$ and by $\Delta \hat{b}$ and $\Delta {\hat{b}}_0$ , respectively. For the initial condition, $\Delta {\hat{\theta }}_0^{-}$ was set as 0 and $\Delta {\hat{b}}_0^{-}$ was set as $\Delta \hat{b}$ at the final measurement. The Kalman filter was applied repeatedly until its outputs converged.

Ant colony algorithm to obtain optimal parameters

In the natural world, ants of some species wander randomly, and upon finding food return to their colony while laying down pheromone trails. If other ants find such a path, they are likely not to keep traveling at random, but instead follow the trail, returning and reinforcing it if they eventually find food. Pheromone evaporation also has the advantage of avoiding the convergence to a locally optimal solution. If there were no evaporation at all, the paths chosen by the first ants would tend to be excessively attractive to the following ones. The overall result is that when one ant finds a good/short path from the colony to a food source, other ants are more likely to follow that path, and positive feedback eventually leads to all the ants following a single path. The idea of the ant colony algorithm is to mimic this behavior with “simulated ants” walking around the graph representing the problem to solve.

The two uncertain parameters of the Kalman filter, for which different values lead to different results, are the angular velocity weight $Q\_\mathrm{gyro}$ and acceleration weight $Q\_\mathrm{acc}$ , and the ant colony algorithm is employed to find their optimal parameter combination in two-dimensional (2D) space. The search criteria of the ant colony algorithm should first be defined. In the proposed method, the search criteria seek to minimize the bias between the reference angles obtained from the standard subject experiment and the angles calculated by processing the IMU data obtained from the human subject experiment.

First, the solution of the optimization problem is determined as $\mathrm{\omega }=({\mathrm{\omega }}_1,{\mathrm{\omega }}_2{)}^{\mathrm{T}}$ and the range of each parameter is confirmed as ${\mathrm{\omega }}_{\mathrm{il}}\le {\mathrm{\omega }}_i\le {\mathrm{\omega }}_{\mathrm{iu}}(i=1,2)$ , where ${\mathrm{\omega }}_{\mathrm{iu}}$ and ${\mathrm{\omega }}_{\mathrm{il}}$ represent the upper and lower bound of each parameter, respectively. Then, each parameter is divided into N intervals, each with a value denoted h, as given by equation (19), to construct a two column by ( $N+1$ )-row matrix

$$h_i=\frac{{\mathrm{\omega }}_{\mathrm{iu}}-{\mathrm{\omega }}_{\mathrm{il}}}{N}$$

(19)

In this model, each figurative ant constructs a solution by choosing a point at each row according to the pheromone value $\tau$ at each point, as depicted in Figure 3.

Figure 3.

Schematic of parameter space meshing.

The probability $P_{\mathrm{ij}}$ of choosing a point at row i and column j, and the initial values of the pheromones, are described as follows

$$P_{\mathrm{ij}}=\frac{{\tau }_{\mathrm{ij}}(t)}{{\tau }_{i0}(t)+{\tau }_{i1}(t)+···+{\tau }_{\mathrm{iN}}(t)}$$

(20)

$${\tau }_{\mathrm{ij}}(0)=\frac{1}{N+1}$$

(21)

where for point $i$ , $j$ , ${\tau }_{\mathrm{ij}}(t)$ is the pheromone value and ${\tau }_{\mathrm{ij}}(0)$ the initial pheromone value under the assumption that the sum of pheromones for each row is 1. Then, equation (21) can be simplified as

$$P_{\mathrm{ij}}={\tau }_{\mathrm{ij}}(t)$$

(22)

After the $m$ ants construct $m$ solutions, these solutions are used as inputs of the objective function, which is the bias between the reference angles obtained from the standard subject experiment and the angles calculated by processing the IMU data obtained from the human subject experiment. The current optimal solution $S_{\mathrm{ib}}$ is obtained from the output of the objective function at each search cycle. When the first cycle is completed, the first $S_{\mathrm{ib}}$ is assigned as the global optimal solution $S_{\mathrm{gb}}$ . The value of $S_{\mathrm{ib}}$ obtained from each subsequent search cycle is then compared with the existing $S_{\mathrm{gb}}$ . If $S_{\mathrm{ib}}$ is better than $S_{\mathrm{gb}}$ , $S_{\mathrm{gb}}$ is replaced with $S_{\mathrm{ib}}$ . Then, the pheromones are updated based on the position of $S_{\mathrm{gb}}$ as follows

$${\tau }_{\mathrm{ij}}(t+m)=(1-\rho ){\tau }_{\mathrm{ij}}(t)+\Delta {\tau }_{\mathrm{ij}}(t+m)$$

(23)

$$\Delta {\tau }_{\mathrm{ij}}(t+m)=\rho$$

(24)

where $\rho$ is the pheromone volatility, and $\Delta {\tau }_{\mathrm{ij}}(t+m)$ is the increased value of the pheromone. The principle of updating is that the pheromones at all points will volatilize, but only the pheromones of the points included by $S_{\mathrm{gb}}$ can be increased.

When the number of search cycles reaches an upper bound $N_{c\max }$ , the points at each row with the largest pheromone values and the number of columns $(c_1,c_2,\dots ,c_n)$ associated with them are obtained. Then, the bounds of each parameter are narrowed as follows

$${\mathrm{\omega }}_{\mathrm{il}}=\{\begin{array}{cc}\hfill {\mathrm{\omega }}_{\mathrm{il}}+(C_i-\Delta )h_i,\hfill & \hfill {\mathrm{\omega }}_{\mathrm{il}}+(C_i-\Delta )h_i\ge {\mathrm{\omega }}_{i0}\hfill \\ \multicolumn{1}{c}{{\mathrm{\omega }}_{i0},}& \hfill \mathrm{else}\hfill \end{array}$$

(25)

$${\mathrm{\omega }}_{\mathrm{iu}}=\{\begin{array}{cc}\hfill {\mathrm{\omega }}_{\mathrm{il}}+(C_i-\Delta )h_i,\hfill & \hfill {\mathrm{\omega }}_{\mathrm{il}}+(C_i-\Delta )h_i\ge {\mathrm{\omega }}_{\mathrm{iN}}\hfill \\ \multicolumn{1}{c}{{\mathrm{\omega }}_{\mathrm{iN}},}& \hfill \mathrm{else}\hfill \end{array}$$

(26)

where $\Delta$ is the set value by which to narrow the bounds, and ${\mathrm{\omega }}_{i0}$ and ${\mathrm{\omega }}_{\mathrm{iN}}$ are the maximum and minimum values of each row, respectively.

After the ranges of each parameter have been updated, the parameter matrix is reconstructed according to equation (19), and the initial pheromone values of each point are reset according to equation (21). The search process is then repeated until the bias between ${\mathrm{\omega }}_{\mathrm{il}}$ and ${\mathrm{\omega }}_{\mathrm{iu}}$ finally achieves an established accuracy $\epsilon$ , as illustrated by the overall process flowchart shown in Figure 4.

Figure 4.

Flowchart of ant colony algorithm.

The optimal parameter ${\mathrm{\omega }}_i^{*}$ is obtained as follows

$$|{\mathrm{\omega }}_{\mathrm{il}}-{\mathrm{\omega }}_{\mathrm{iu}}|<\epsilon$$

(27)

$${\mathrm{\omega }}_i^{*}=\frac{{\mathrm{\omega }}_{\mathrm{il}}+{\mathrm{\omega }}_{\mathrm{iu}}}{2}$$

(28)

The computational complexity of the proposed algorithm can be analyzed according to the above algorithm steps. The steps involving the highest computational complexity are given as follows: (1) the structure of the solution involving the calculation of the objective function value, (2) the process of updating pheromones, and (3) the space subdivision. The first step represents the highest computational complexity and can therefore represent the computational complexity of the overall algorithm, which is given as follows

$$O(s·m·(N+1)·n·N_{c\max })$$

(29)

where $s$ represents the number of the search cycles when the precision reaches the requirement, and $n$ represents the number of parameters, which is 2 here.

Results

Data acquisition

The data-acquisition hardware was employed to collect signals from a standard subject and a human subject. The standard subject was a servo motor, and the human subject was a healthy 20-year-old male. Data collected from the standard subject were treated as the original signal. Data collected from the human subject were treated as the result of combining the signal and noise. The Kalman filter was used to fuse the acceleration and angular velocity data to determine the angle. Integration of the angular velocity represented the change in the angle. Prior to data acquisition, all IMUs were calibrated to minimize the offset. However, even a precise calibration yields a degree of non-zero bias, which inevitably imparts a degree of drift to the angles integrated from the gyroscope data. In this study, the acceleration was used to rectify the steady-state error. The detection noise of the accelerometers was mainly derived from vibration, leading to a velocity detection bias. The accelerometer noise was defined as the resultant acceleration obtained from human subject measurements subtracted by the acceleration due to gravity. The signals from the standard subject were treated as the ground truth. Then, the noise obtained from the human subject was added to the ground truth data to simulate a realistic standard signal. A quaternion coordinate system, widely used in the inertial navigation field, was employed to deduce the rotation angles based on the angular velocity data. As ideal control parameters, the coordinate system represents a rigid body’s Euler axis and rotation angle during the attitude transformation. We subsequently obtained the rotation angles based on the acceleration data by calculating the angles between the segment vector in the local sensor coordinate system and the gravity vector. After obtaining these two types of rotation angles, the two angles were fused by means of the Kalman filter. Ant colony algorithm optimization was employed to obtain an optimal combination of two of the parameters required by the Kalman filter to achieve an optimal result relative to the ground truth. Experiments were conducted for both the standard subject and the human subject, and the results demonstrated the feasibility of the proposed algorithm.

Two different experiments were designed. In the first, data were acquired from the standard subject to obtain a series of reference rotation angles for implementing the ant colony algorithm. Data were then acquired from the human subject to verify the feasibility of the algorithm. In both of these two experiments, the acquired IMU data were stored in the SD card, which was convenient for offline processing.

In the standard subject experiment, a motorized rotation stage was placed on an optical platform, keeping the rotation plane vertical to the horizontal plane. A single IMU sensor was mounted on the rotation plane. At the beginning of data acquisition, the servo motor is held stationary for 2 s and then rotated clockwise 2° over a period of 0.17 s, held stationary for 0.73 s and then rotated clockwise 90° at a constant angular velocity over a period of 2.7 s, held stationary for 0.6 s and then rotated clockwise 2° over a period of 0.17 s, held stationary for 0.73 s and then rotated anti-clockwise 94° at a constant angular velocity over a period of 3.56 s, and held stationary for 3 s before the end of data acquisition. Over the course of the experiment, the acceleration and angular velocity data obtained from the IMU were collected, and the rotation angles were recorded as the ground truth of the experiment. The experimental record of the rotation angles is given in Figure 5.

Figure 5.

Record of programmed motor rotation angles.

During the human subject experiment, the IMU was mounted on the subject’s thigh segment, and the data-acquisition electrical circuit board was bound at the subject’s rear waist. The subject was first instructed to maintain an upright posture for 5 s, keeping the lower limb segment vertical for the sensor to acquire the direction of gravity for calibration purposes. The subject was then instructed to walk on level ground at a medium speed while IMU data were sampled. The resultant acceleration subtracted by the acceleration due to gravity was treated as the acceleration noise. After collecting the acceleration and angular velocity data of the IMU, the accelerometer noise was added to the standard subject signal to simulate a realistic standard signal.

Verifying algorithm performance

After the reference rotation angles were obtained from the standard subject experiment, the ant colony algorithm was applied to improve the performance of the sensor fusion by searching the optimal parameters in the Kalman filter. To verify the performance of the ant colony algorithm in terms of the searching efficiency or the optimality of the search results, we compared its performance with that of two different algorithms.

Comparison with traversal method

In the Kalman filter there are two uncertain parameters, the angular velocity weight $Q\_\mathrm{gyro}$ and the acceleration weight $Q\_\mathrm{acc}$ , which we use in the ant colony algorithm to find their optimal parameter combination in 2D space. Some other methods can also be used to find the optimal parameter combination, such as the traversal method. Depending on the data acquired from the standard subject experiment, as shown in Figure 2, we compared the operational cycling of the ant colony algorithm with that of the traversal method. For example, in the test experiments, we restricted the parameters as $Q\_\mathrm{gyro}\in [0,0.8]$ and $Q\_\mathrm{acc}\in [0,0.08]$ . To achieve an equivalent precision, such as 0.001, the operational cycling of the ant colony algorithm required only 3761 iterations on average, but the traversal method required over 60,000 iterations, and the computing time of the ant colony algorithm was only one-third that of the traversal method. According to the operational cycling and computing times of these two algorithms, the ant colony algorithm demonstrated a higher search efficiency than the traversal method.

Comparison with use of empirical parameter values

In most Kalman filter applications, the parameters are assigned empirical values. In the present circumstance, most researchers assign the parameters as $Q\_\mathrm{gyro}=0.03$ and $Q\_\mathrm{acc}=0.01$ . For an adequate comparison, the sum of the squared errors for our method, ${\psi }_{\_a}$ , and that of the empirical method, ${\psi }_{\_e}$ , were defined as

$${\psi }_{\_a}=\sum _{i=1}^n{e}_i^2,e_i={\theta }_i-{\theta }_i^a$$

(30)

$${\psi }_{\_e}=\sum _{i=1}^n{e}_i^2,e_i={\theta }_i-{\theta }_i^e$$

(31)

where ${\theta }_i$ is the dip angle calculated by a video system, which was used as the benchmark angle, and ${\theta }_i^a$ and ${\theta }_i^e$ are the dip angles calculated by our algorithm and the empirical method that did not find the optimal parameters, respectively. The values of ${\psi }_{\_a}=1.9305{\mathrm{rad}}^2$ and ${\psi }_{\_e}=6.7723{\mathrm{rad}}^2$ were obtained from test experiments. According to the sum of squared errors, the ant colony algorithm provided much more accurate results than the empirical method.

Human experiment based on proposed algorithm

The optimal values of $Q\_\mathrm{gyro}$ and $Q\_\mathrm{acc}$ obtained by the ant colony algorithm were employed in a human experiment involving a basketball player from Xiamen University, China, as the experimental subject. The experimental subject was 21 years old, played the center position, was a dominant player, and was 2.03 m in height and 100 kg in weight, with a thigh length of 58 cm and a shank length of 50 cm. The height of the subject’s foot was 11.5 cm and its length 32 cm. IMU sensors were mounted on the subject’s body, and a single IMU was mounted on the subject’s left thigh. The data-acquisition electrical circuit board was bound at the subject’s rear waist. The initial posture maintained the lower limbs in a vertical orientation to calibrate the sensor with respect to the direction of gravity. Subsequent motion involved dribbling a basketball with the subject’s right hand (i.e. the predominant hand) while walking in straight line. Video was captured simultaneously with a sampling frequency of 24 Hz. Figure 6 shows a series of continuous frames obtained over a single gait cycle period. Figure 7 shows the deduced results, where the horizontal axis represents the time, and the vertical axis represents the angle between the thigh segment and the gravitational direction. A comparison of the calculated results with the video verified the feasibility of the proposed system, including the hardware and the algorithms employed.

Figure 6.

Continuous video frame series representing a single gait cycle period.

Figure 7.

Deduced left thigh dip angle over a single gait cycle period.

Discussions and conclusion

In this paper, we presented a lower limb segment angle measuring system employing a six-axis IMU with a 3D accelerator and 3D angular velocity outputs. To eliminate gyroscope drift, the Kalman filter was utilized to fuse the respective segment dip angles deduced from acceleration and angular velocity measurements. The Kalman filter employs two uncertain parameters that are typically assigned empirical values. In the proposed method, measurement results obtained from standard subject and human subject experiments were employed by the ant colony algorithm to find the optimal parameter value combination that led to a minimum bias from the ground truth. The human experiment involved a basketball player as the test subject, where the dip angle of the subject’s thigh during the experiment was calculated offline. The segment angles were then deduced from the gravitational direction. The experimental results demonstrated that the wearable sensors employed can provide quantitative measurements of human gait motion, such as joint angles. The accuracy of an optical tracking system is greatly dependent on the position of the installed optical markers, and it is very difficult to eliminate bias. ³³ However, in the servo motor reference system, the programmed rotation angles are highly accurate.

The proposed lower limb segment angle measuring system still requires improvement, and this represents ongoing future work. Because the method employs no magnetic field sensing, the yaw angles cannot be obtained, which causes difficulties when calculating the joint angle. To solve this problem, the angle reversion system is being improved. A nine-axis IMU sensor with magnetic detection will be employed to upgrade the system, and orientation reversion algorithms other than the Kalman filter will be implemented to test their respective advantages. Ultimately, a real-time lower limb angle measurement system will be developed for rehabilitation engineering product applications.