MEMS Inertial Sensors Based Gait Analysis for Rehabilitation Assessment via Multi-Sensor Fusion

Gait and posture are regular activities which are fully controlled by the sensorimotor cortex. In this study, fluctuations of joint angle and asymmetry of foot elevation in human walking stride records are analyzed to assess gait in healthy adults and patients affected with gait disorders. This paper aims to build a low-cost, intelligent and lightweight wearable gait analysis platform based on the emerging body sensor networks, which can be used for rehabilitation assessment of patients with gait impairments. A calibration method for accelerometer and magnetometer was proposed to deal with ubiquitous orthoronal error and magnetic disturbance. Proportional integral controller based complementary filter and error correction of gait parameters have been defined with a multi-sensor data fusion algorithm. The purpose of the current work is to investigate the effectiveness of obtained gait data in differentiating healthy subjects and patients with gait impairments. Preliminary clinical gait experiments results showed that the proposed system can be effective in auxiliary diagnosis and rehabilitation plan formulation compared to existing methods, which indicated that the proposed method has great potential as an auxiliary for medical rehabilitation assessment.


Introduction
Human walking contains important physiology, kinematic and dynamic information. There are many application prospects of human gait analysis in real life, such as monitoring the patient's recovery progress in clinical practice, the control strategy of bionic robots, etc. For instance, stroke is a common neurodegenerative condition with a principal symptom of progressive limbs movement disorder. The prevalence of this neurological disease has been estimated at 120 affected individuals per 100,000 and 75% of the victims suffer from sequelae (mostly involve gait dysfunction) afterwards according to epidemiological statistics [1][2][3]. An insightful example of gait disorder can be observed in stroke-note that this study is not merely focused on stroke. However, a great majority of subjects who took part in the gait analysis experiment suffer from stroke. Furthermore, stroke patients share most gait symptoms with common neurological disorders in clinical practice. Specifically, stroke is a chronic neurological disorder associated with hemiplegia, lack of balance and abnormal gait. Therefore, objective and accurate inspection of gait parameters (as shown in Table 1) is of great help to a neurologist for appropriate assessment and diagnosis of stroke patients. To this end, the findings of this study could be helpful for revealing the pathogenesis of gait disorders. Stroke normally results in gait asymmetry, which can be reflected in the differences of gait phase partition between two feet due to paresis and imbalance. To sum up, it is crucial to discover the main components of gait disorders.
In view of this situation, it is critical to develop an objective and quantitative approach to assess the patients' physical condition. This paper established a wearable gait analysis platform based on a MEMSsensor and body sensor network. The platform can be used to collect acceleration, angular velocity and the geomagnetic signals in the process of walking movement. Accurate gait parameters can be calculated through a sensor data fusion algorithm and error correction process. Table 1. Typical spatio-temporal gait parameters.

Gait Parameter Description
Stride length (m) Distance between two consecutive footprint of the same foot.
Stride speed (m/s) Stride length divided by walking cycle.
Stride frequency Number of steps taken per minute during walking.
Walking cycle (s) Duration of a single stride, inversely proportional to cadence.
Stance time (s) Duration of stance phase when feet contact with the ground, starting with initial-contact (IC) and ending with foot-off (FO) of the same foot.
Swing time (s) Duration of swing phase when feet swing above the ground, starting with FO and ending with IC.
Clearance (m) Foot elevation in swing phase, which reflects the muscular strength of lower limbs and can be diversified as maximum and minimum foot elevation.
Plantar & dorsiflex (degrees) The angle between the dorsum of the foot and the back of the leg.
Knee ROM (degrees) Range of knee flexion during a single stride.
With the rapid development of modern medical technology, the concepts of medical treatment have been gradually changed to "Prevention first". In this case, it is quite necessary to conduct the acquisition and processing of health information in advance, so that early medical diagnosis and intervention are feasible. Meanwhile, in order to implement ambulatory monitoring without affecting the subjects' normal physiological activities, the traditional wire communication mode gradually shifts to a wireless mode, which is likely to be more persuasive. Moreover, innovative sensing technology is indispensable for continuous monitoring since the medical monitoring equipment has been developed towards microscale and a long-term span. To this end, the emerging Body Sensor Network (BSN) might serve as a peripheral node of the Internet of Medical Things or even the ubiquitous network.

Related Works
Quantitative gait analysis system mainly including a camera system [4][5][6], electromyography measuring system [7] and force platform [8,9]. The camera system consists of multiple high resolution cameras located in an indoor space, the orientation and position information of the target subject can be calculated using attached highlighting reflective spots. The electromyography measuring system detects human lower limb muscle signals by surface electromyography in the waking process; force platform reflects the change of plantar pressure during walking. However, the applications of above gait analysis system are limited in clinical practice, and the main reasons lie in three aspects. Firstly, the systems are expensive, which might be a barrier to routine use. Secondly, the usage of the systems are complex and require special operation, and it usually takes hours to complete the whole gait measurement process. Finally, specific space is normally needed to perform gait analysis using the above systems. In particular, the camera system may need more than 100 hundred square meters [10][11][12][13]. Table 2 lists a brief comparison of mainstream gait analysis methods.
Considerable research has been conducted into the progression of gait dysfunction through the various stages of stroke. Specifically, stroke subjects experience decreased stride length, cadence and walking speed, significant variability in stride length and gait cycle, and Walking imbalances [10,[14][15][16]. Chang et al. [17] employed a specialized wearable system and found that stroke subjects demonstrated decreased gait velocity, stride length and prolonged double support phase. They further identified a high correlation between these gait parameters and age of onset. Further investigation into gait impairments in subjects with neurological diseases has also indicated that the degree of gait abnormality and the disease progression. Previous research has highlighted the advantages of quantitative gait analysis in gait diagnosis; however, laboratory-based systems such as optical tracking and plantar pressure measurement are typically expensive and are not available in ordinary clinical settings [18][19][20]. Therefore, significant interests have increased rapidly in the development of alternative gait analysis tools. With the maturity of microelectromechanical systems and the development of information fusion technologies, the application of inertial motion analysis technology is becoming more and more extensive [8,[21][22][23][24][25]. Due to the noticeable advantages of small size and low cost, wearable sensors can be mounted directly on the body segment with no need for specified test environment [24,[26][27][28]. Such system may also serve as a good supplement of the gold standard including optical system and plantar pressure monitoring system. In previous studies, we have adopted a wearable inertial sensor in a walking distance calculation and walking pattern classification [29][30][31]. Ambulatory measurement of the participant's trunk inclination using inertial measurement unit (IMU) was carried out by Farris et al. [3]. Bao et al. [32] developed a smart shoe for gait analysis using force sensitive resistors and IMU sensors. Luinge et al. [33] proposed the estimation of arm orientation by wearable inertial sensors. Dejnabadi et al. [34] introduced an approach to accurate measurement of joint angles based on IMU. However, due to their inability to detect heading reference, inertial based systems generally fail to measure differential orientation, a prerequisite for computing the 3D knee flexion angle recommended by the Internal Society of Biomechanics [35]. Roetenberg et al. [36] developed an ambulatory position and orientation tracking method fusing magnetic and inertial sensing. Since magnetometers measure the strength and direction of the local magnetic field, the geographical north direction can be found. In this case, the initial heading orientation can be obtained with the supplement of a magnetometer. Moreover, the system remains self-contained, which means it does not rely on any external infrastructure [37]. In addition, there are already wireless IMU BSN commercial products such as Trigno Avanti (Delsys Inc., Natick, MA, USA), Mvn Suit (XSens Inc., Enschede, The Netherlands), Perception Neuron (Noitom Inc., Beijing, China) and iSen (STT Systems Inc., San Sebastian, Spain). The current limitations of the state-of-the-art mentioned in the literature are the sensor alignment and integral error. Cost-effectiveness and stability are two other concerns. Moreover, little research about the follow-up monitoring of patients' lower limbs has been carried out. Therefore, the contributions of this paper include the sensor alignment method and the availability of follow-up monitoring of patients' key gait parameters.
The rest of this paper is organized as follows: Section 3 describes the structure of the proposed gait analysis system and the methodology used to estimate the gait parameters during walking; experimental results are given in Section 4; and the potential applications of gait analysis are discussed in Section 5, which concludes the paper as well.

System Setup
An ambulatory gait analysis system has been developed based on IMU. We have assembled the IMU from accelerometers/gyroscopes chipsets including ST Microelectronics (Geneva, Switzerland) and Analog Devices (Boston, MA, USA). A multi-sensor fusion algorithm is used to estimate gait parameters. The inertial measurement unit (IMU) can be strapped on both feet, shank and thigh via adjustable elastic straps with hook & loop. We chose straps over housing shells due to their flexible structure, strong adaptability, lower cost and good maneuverability. The installation is quite simple, which can be finished in several minutes. We have designed a fastener with which the sensor nodes can be fixed on the specified location firmly, ensuring the estimation accuracy. No special laboratory is needed, and the gait measurement can be performed just in the corridor or in the ward. The patients can even keep follow-up monitoring in the community after they discharged from the hospital. The performance-to-price ratio is relatively high and it is convenient to automatically generate a gait diagnostic report. Moreover, useful contrastive analysis can be made with repeated gait analysis. The principle and structure of the proposed gait analysis system is shown in Figure 1.
The sensor array includes triaxial accelerometers, gyroscope and magnetometer. Sensor performance specification is shown in Table 3. Raw sensor data is logged at 100 Hz and then forwarded to a receiver via a 2.4 GHz wireless network. The actual linear motion acceleration is used to calculate position. Note that the gravity component (g = [0, 0, 9.81] T m/s 2 ) is normally eliminated from the resultant acceleration signal before estimating the position by integral operation of actual motion acceleration two times.

Accelerometer Non-Orthogonal Error Estimation
In general, the accelerometer orthogonal error is less than 1 • , which is normally indicated in the product technical manual. Orthogonal errors exist in three-axis, as shown in Figure 2, and the non-orthogonal error can be expressed as follows: Formula (1) can be simplified using the approximation of the trigonometric function value and then we have the following equation: The accelerometer merely senses gravity under stationary state. A relationship exists between accelerometer observations g s and the true gravitational acceleration ±bg as follows: Meanwhile, g T s E T Eg s g s 2 = 1. After establishing the non-orthogonal error correction model, the ellipsoid fitting method is introduced to calculate the non-orthogonal error angle. We can get three components of gravity vector ±bg with respect to three sensitive axes as follows: while the accelerometer observations G s = [a x , a y , a z ]T should satisfy the following equation theoretically: which means the distribution of measured values is a spherical with the radius of g g g 2 ; however, actually, the accelerometer measurement distribution is an ellipsoid due to the existence of nonorthogonal errors. Set O = [a, b, c, f , g, h, p, q, r, d] T and quadric equation can be written as follows: Then, we can get the coefficient matrix from Equation (4) after ellipsoid fitting, Define Then, Equation (7) is the description of ellipsoid surface. Select n group of accelerometer measurements, [x i , y i , z i ], i = 1, · · · , n, denote C = [X 1 , X 2 , · · · , X n ] and ] T , and then the ellipsoid fitting is converted into the following constraints: Define the coefficient matrix as follows: Lagrangian function can be introduced to convert the constraint problem to the following equation: It has been proved that this constraint problem has a unique solution in the field of mathematics [38]. Therefore, the ellipsoid coefficient vector can be determined, and then the ellipsoid radius and the nonorthogonal error angle can be calculated based on the ellipsoid coefficient. Define the symmetrical coefficient matrix and the translation vector S T = [2p, 2q, 2r] T , Then, the transformation H = N T H + R can convert the quadric equation into a standard ellipsoid equation: x 2 a 2 + where N is the eigenvector of matrix S, denote D as the main diagonal matrix of S, and then we have S = NDN T and R = −(2D) −1 NT. Ellipsoid radius can be obtained by the following formula: After calculating the ellipsoid coefficient, we can calculate the nonorthogonal error angle by the analytic method. Equation (3) can be used to compensate the non-orthogonal error of the accelerometer.

Stance Phase Detection by Decision Level Data Fusion
Sensor drift is an inherent property that results in linear growing integration errors in attitude and position estimation. In particular, position errors grow proportional to the square of the acceleration error. To this end, the widely used Zero Velocity Updating (ZVU) method is adopted in this paper. Though the valid interval of ZVU algorithm is illustrated in literature [31,[39][40][41], the effectiveness of the ZVU technique largely relies on the stance phase detection as shown in Figure 3. This paper takes two criteria to determine stance phase in each gait cycle. These key periods are determined by calculating the squared Euclidean norm of acceleration values, as shown in formula (16): where a x , a y and a z represent the triaxial acceleration measurements of foot sensors: whereS N is the mean of S i over N samples. Meanwhile, angular rate energy E gyro is adopted as the other criterion (18). The second moment is used to detect stance phase, which is defined in the following formula [39]: where W is the window size selected according to the sensors sampling rate; ω i = [ω x,i , ω y,i , ω z,i ] T is the triaxial angular velocity vector; and σ 2 ω is the gyroscope noise variance. λ 1 and λ 2 are empirically predefined thresholds. The detection results ( R) are sequences consisting of "zero" and "one". The algorithm continually finds the interval when ZVU is valid and updates the corresponding v G (t) as [0, 0, 0] T based on the two indicators above.

Knee Angle Estimation
Knee flexion angle can be calculated by fusing the multi-sensor data. Figure 4a demonstrates the calculation principle and Figure 4b,c show the calculated swing angles and knee flexion based on the gyroscope observations. In addition, the calculation method of knee flexion angle is as follows: where θ 0 thigh and θ 0 shank are initial angles between lower limbs (thigh and shank, respectively) and the gravity direction when the subjects are standing still at ease, which can be calculated by the measurements of the accelerometer as proposed in the previous study [29,30]. ω thigh and ω shank are angular velocity values of thigh fixed and shank fixed sensors.

Attitude Estimation and Quaternion Correction in IMUs via Sensor Fusion
IMUs refer to sensor modules consisting of three-dimensional accelerometers, gyroscopes and magnetometers (in some cases magnetometers are not included). According to physical and dynamical theory, acceleration measurements can be integrated once to acquire linear velocity and twice to obtain relative position change based on the previous observations. Likewise, angular velocity observations from gyroscopes can be integrated once to estimate the attitude change between two consecutive measurements. In practice, inertial sensors are prone to be disturbed by system noises, drift and measurement errors. All these factors would cause significant integration errors when the raw sensor data are used for integration. In some cases, when the subject moves slowly or stays static, one can merely adopt accelerometers or inclinometers to directly determine the 3D attitude with acceptable results, which avoids the integration operation. However, it is inevitable that external acceleration applied to the accelerometers could ruin the attitude detection based on the calculation of gravity vector components in these orthogonal planes. In most cases, multiple sensor data fusion is necessary to determine 3D attitude. Since the collected data are discrete, we need to perform interpolation between q m and q n and the interpolation principle is shown in Figure 5.
Simple linear interpolation is valid in some cases. However, it can not effectively describe the curve between q m and q n . To ensure that the angle θ between q m and q t is linear, i.e., θ(t) = (1 − t)θ + tθ, we then choose slerp function slerp(q m , q n , t) to conduct smooth interpolation of quaternions as follows: Normally standardized operation is necessary: Moreover, the quaternion number can also avoid the singular point problem of Euler angle representation [27,42], when the pitch angle reaches ±90 • . The three-dimensional attitude represented by quaternions are: The attitude of the updated sensor can be obtained by solving the differential equation of quaternions. The differential equation of quaternions can be expressed as: On the basis of the widely used sensor coordinate transformation, we can get the updated attitude quaternion with inevitable errors. As for a certain vector, its magnitude should be the same though it is expressed in different coordinate frames. The magnitude deviation caused by coordinate transformation can be adopted to adjust the rotation matrix. This paper uses two reference vectors (gravity vector and magnetic vector) to modify the quaternion. In the static state (without linear acceleration), gravity vector [0, 0, 1] T is converted into [c x , c y , c z ] T after coordinate transformation, while accelerometer observations are [a x , a y , a z ] T . Then, [c x , c y , c z ] T and [a x , a y , a z ] T both represent the gravity vector in the sensor frame. Then, we can obtain error matrix e s g by multiplying both vectors: The error matrix can be used to correct the attitude quaternion. Proportional-Integral (PI) feedback adjustment is introduced hereby, where ω t is a three-axis angular velocity component, which can be used to correct quaternions combined with skew symmetric matrix. k p and k i are proportion coefficient and integral coefficient, respectively. Both parameters are ascertained after repeated experiments on different specific groups. Moreover, in practice, the measurement error of MEMS magnetometer is non-negligible, the magnetometer measurement error mainly includes environmental interference and inherent error [12,37]. Researchers have conducted various calibration methods using high precision instruments and equipments [38,43,44], however, these methods are time-consuming and the calibration effect largely relies on the precision of equipment. This work uses a calibration method which allows non-experts to easily implement the calibration procedure. To calibrate the magnetometer, sensor nodes were rotated along "Figure of eight knot" trajectory for several times before they were mounted on human lower limbs. Then the outputs were wirelessly collected and be used to perform the magnetometer correction. The correction of magnetometer should be performed near the body segment where the sensor was mounted in an indoor experimental scene [45]. In conclusion, the flowchart of the proposed gait parameters estimation approach can be summarized in Figure 6.

Experimental Results
The proposed gait analysis system has been adopted to carry out preliminary clinical trials at the Second Affiliated Hospital of Dalian Medical University. Gait data used in this study consist of walking stride time series from 30 healthy adults (22-45 years old), 20 patients with stroke (46-77 years old), and 20 patients with joint disease (30-58 years old). All subjects were instructed to walk continuously on level ground along an obstacle-free corridor for more than 15 m. Note that we have conducted walking trials including U-turns and stair climbing into the evaluation and the proposed method works well. Straight walking is a simplified case suggested by the clinician, and straight walking is mostly adopted in the observational method in clinical practice.
According to the obtained gait parameters in Table 4, the stride lengths, stride speed and feet clearance are relatively low in both neurological and arthropathy patients, which are also consistent with clinical observation. Similarly, mean values of stance time associated with patients are significantly higher than those of healthy subjects. The results from statistical results presented in the table indicate strong evidence of the capability of the typical gait parameters in characterizing the walk of healthy subjects and patients. In this regard, these gait-related symptoms can be explained by clinicians for diagnostic and treatment purposes because the disease progression can be quantificationally monitored via these observations.  Table 5 presents the maximum of joint angle during walking. The knee flexion angle is constantly positive; the positive ankle joint angle represents dorsiflexion and a negative value signifies plantarflexion. In the course of stance phase, the knee flexion angle increases; meanwhile, the ankle dorsiflexion turns to plantarflexion; note that the maximum ankle plantarflexion appears in the final stage of stance phase; in the course of swing phase, knee flexion angles reach the maximum while ankle plantarflexion turns to dorsiflexion, followed by the next stance phase. Knee flexion range of motion (ROM) is often evaluated using a goniometer in rehabilitation clinics or in hospital wards. The more knee ROM regained during the therapeutic process, the better knee recovery would be affirmed and the sooner early discharge could be guaranteed. In this study, we conducted data collection and knee ROM analysis on 30 healthy subjects, 20 neurological patients (mainly stroke patients) and 20 arthropathy patients, respectively. The first data collection is performed before medical treatments, while the remaining data collections occurred at two weeks and six weeks after treatments, respectively. Figure 7 illustrates the knee ROM recovery status of one typical arthropathy patient who received minimally invasive surgery and one typical stroke patient undergoing rehabilitation training, respectively. The results show that both patients recovered significantly in terms of knee ROM after receiving six weeks of medical treatments. By six weeks after minimally invasive surgery, the knee ROM of arthropathy patient almost returns to the normal range and the gait symmetry is much better when the pains were alleviated.
The one-way analysis of variance (ANOVA) results of knee ROM between bilateral lower limbs are shown in Tables 6 and 7. We can conclude that patients showed large standard deviations in knee ROM, which is a significant feature different from healthy subjects. In this paper, the hypothesis is the symmetrical (balanced) bilateral knee angle ROM, p-value < 0.05 should be interpreted as the hypothesis is true, and the hypothesis is invalid for subjects with p-value > 0.05. Results showed that, as for arthropathy patients, no significant knee ROM difference exists on bilateral lower limbs after six weeks' treatment based on the ANOVA analysis results (p-value = 0.0046); however, the stroke patient still has significant knee ROM different on bilateral lower limbs after six weeks of treatment (p-value = 0.8637). In fact, many stroke patients still have obvious asymmetry between bilateral knee ROM even after several months, though the symptoms may be relieved to a great extent.

Feet Clearance Monitoring
According to Figure 8a, it is observed that there exists foot elevation asymmetry (left foot: ∼15 cm, right foot: ∼5 cm) of a stroke patient with hemiplegia symptoms on the right foot, which is consistent with clinical observation, indicating that the patient loses the ability to keep balance. Figure 8b presents the foot elevation asymmetry of an arthropathy patient, which presents the gait walking disorder from the perspective of feet elevation statistics.

Discussion and Conclusions
This paper aims to provide ambulatory and robust measurements of human gait, and we adopted body-worn sensors to estimate gait parameters. Digitalized and objective gait information can act as desirable guidance for making and adjusting rehabilitation plans, and the results of the preliminary clinical gait analysis experiments have also verified the accuracy of this method for human limbs motion capturing. With no pressure sensor for the stance phase detection and no optical device for integral error elimination, the pure ZVU-aided gait analysis system using body-worn IMU can achieve a good auxiliary diagnosis performance.
Experimental studies have been presented for an Magnetic Angular Rate and Gravity (MARG) unit with reference measurements obtained via a precision optical measurement system, i.e., the NDI Polaris Spectra System (Northern Digital Inc., Waterloo, ON, Canada). The proposed gait analysis system accuracy was validated and the three-dimensional position estimation error is less than 0.015 m, as shown in Figure 9. A comparison experiment with an optical system has demonstrated the accuracy and feasibility of the proposed principle of error correction, and the results of the preliminary clinical gait analysis experiments have also verified the accuracy of this method for human limbs motion capturing.
The current work provides new insights to better understand the biomechanics of walking due to neurological diseases. In addition, they appear to be valuable tools that can highlight differences in gait dynamics with respect to stroke patients. In this regard, measurement of gait may possibly afford pertinent clinical information on neuromotor conditions, characterization of some neurological disorders, and rehabilitation. A better understanding of gait differences based on etiology of amputation or fall history may provide useful information to help guide prosthetic prescription or rehabilitation interventions.
(a) (b) Figure 9. System accuracy validation by an NDI Polaris Spectra System (a) sensor placement and reflection points of optical system; (b) error statistics of a three-dimensional thigh sensor position estimation for random lower limbs' movement.