Gait analysis algorithm for lower limb rehabilitation robot applications

When patients with lower limb dyskinesia use robots for rehabilitation training, gait parameters are of great significance for disease diagnosis and rehabilitation evaluation. Gait measurement is usually carried out by using optical motion capture systems, pressure plates and so on. However, it is difficult to apply these systems to lower limb rehabilitation robots due to their high price, limited scope and wearing requirements. At the same time, most of the current applications in robots focus on the basic gait parameters (such as step length and step speed) for robot control or user intention recognition. Therefore, this paper proposes an online gait analysis algorithm for lower limb rehabilitation robots, which uses a lidar sensor as the gait data acquisition sensor. The device is installed on the lower limb rehabilitation robot, which not only avoids the problems of decline in the detection accuracy and failure of leg tracking caused by lidar placement on the ground, but it also calculates seven gait parameters, such as step length, stride length, gait cycle and stance time, with high precision in real time. At the same time, the walking track of the patient may not be straight, and the lidar coordinate system is also changed due to the movement of the lower limb rehabilitation robot when the patient moves forward. In order to overcome this situation, a spatial parameter-splicing algorithm based on a time series is proposed to effectively reduce the error impact on gait spatiotemporal parameters. The experimental results show that the gait analysis algorithm proposed in this paper can measure the gait parameters effectively and accurately. Except for the swing time and double support time, which are calculated with large relative errors due to their small values, the relative errors of the remaining gait parameters are kept below 8 %, meeting the requirements of clinical applications.

1 Introduction

In recent years, lower limb rehabilitation robots have been widely used in rehabilitation training due to their advantages of task repeatability, motion controllability, measurement reliability and high feedback (Qian and Bi, 2015). They mainly helps patients with lower limb dyskinesia from conditions such as stroke, Parkinson's disease, spinal cord injury and so on to carry out dependent ambulation and gradually improve their gait. In addition, the analysis of patients' gait under the lower limb rehabilitation robot can provide objective and accurate evaluation parameters for rehabilitation doctors. The study has shown that quantitative gait analysis is helpful in the early diagnosis of disease and detection of the patients' rehabilitation status in clinical applications. Quantitative data will also enhance their confidence in recovery (Prakash et al., 2018). Traditionally, most hospitals adopt the visual observation method, in which doctors directly observe the walking state of patients and make a subjective judgment of the patient's condition according to their experience. Compared with quantitative analysis based on instruments, this is subjective and has low retest reliability.

With the development of sensor technology, quantitative gait analysis has gradually attracted extensive attention and research. At present, methods of gait feature extraction are mainly divided into several categories: collecting image data of human movement for detection and tracking, feature extraction identification and the establishment of human models (Ahmed et al., 2015; Hynes et al., 2021; Wagner et al., 2023); attaching external markers captured by optical systems (Hutchison et al., 2015; Matjačić et al., 2017; Ziegler et al., 2020) and capturing gait information such as ground reaction force through force plates with dense and high-frequency pressure sensor arrays (Costilla-Reyes et al., 2019; Singh et al., 2017; Vera-Rodriguez et al., 2011). However, the devices used in these methods are often not easily applied to the general public due to their high costs, high hardware requirements, large space requirements and limited scope. Instead, some scholars use wearable sensors of relatively small size to analyze gait, for instance, integrating and attaching inertial sensors (accelerometers, gyroscopes, etc.) to the body surface (thighs, back or feet, etc.) and measuring acceleration and angular velocity changes during walking (Buckley et al., 2020; Gurchiek et al., 2020; Q. Ji et al., 2018; Li et al., 2020); embedding sensors (pressure sensors, accelerometers or gyroscopes) into shoes, socks or insoles (Amitrano et al., 2020; Eizentals et al., 2018; Wang et al., 2019) and integrating sensor systems to construct three-dimensional human joint motion (Baghdadi et al., 2018; Li et al., 2018; Liu et al., 2011). Nevertheless, compared with problems such as poor universality, high wearing requirements and complex operation of wearable sensors, recently, some researchers have put forward the method of measuring gait parameters with lidar sensors as they are low cost and easy to set up and move (Duong and Suh, 2020; Iwai et al., 2019; Yorozu et al., 2015a, b).

Gait analysis has been widely used in assisting robotic control strategies (Cifuentes and Frizera, 2016; Shi et al., 2020; Ye et al., 2020), detecting abnormal gait and assisting machines to ameliorate patients' gait (Lim et al., 2016). Among them, regarding the aspect of gait analysis using lidar, because the sensor itself is characterized by small size and non-contact measurement, it is very suitable for installation on mobile devices. At the moment, there is much research on lidar in mobile robots (Piezzo et al., 2017; Yorozu and Takahashi, 2016, 2020) and intelligent walkers (Bayón et al., 2017; Xu et al., 2018).

In Yorozu and Takahashi (2020), the lidar sensor is mounted on the service robot to track the gait and estimate the pedestrian's body direction, such as forward, backward, U-turn and other movements, namely position, speed and leg state. In Yorozu and Takahashi (2016), the proposed leg detection method can effectively identify the cross steps, reduce the possibility of losing track of legs and better serve the safety collision avoidance of autonomous mobile robots. It realizes leg detection through leg pattern recognition and data association analysis. In Piezzo et al. (2017), a service robot for monitoring and encouraging the elderly to walk is proposed, and the step length of the elderly is measured by lidar. In Xu et al. (2018), a walking-aid robot system is proposed that uses a force-sensing resistor and a laser ranging finder to detect the intended human motion expressed by the user's upper and lower limbs so as to detect if the user is going to fall. A rehabilitation robotic platform for patients with cerebral palsy is proposed in Bayón et al. (2017). Among them, the system is based on the motion control of detecting the user's leg position and motor movement to control the linear speed of the platform, so as to improve the dependence of patients in the rehabilitation process.

In summary, most of the cases applied to robots are focused on gait serving robot control, user intention recognition and gait cycle identification based on offline models. However, the current gait analysis method is rarely applied to lower limb rehabilitation robots. Traditional detection methods also have a variety of problems: due to the complexity of the lower limb rehabilitation robot structure, there is usually an occlusion problem. At the same time, the test subjects need to wear external markers and ensure that the markers will not loosen or fall off easily during the test, which is not suitable for daily walking training of patients. As walking training with rehabilitation robots usually requires long distances, some sensors (such as inertial sensors and floor pressure sensors) cannot meet the measurement requirements well. Based on this, this paper selects a non-contact lidar sensor to collect data, which is installed on the lower limb rehabilitation robot. Without affecting daily rehabilitation training, it can also avoid the problems of decreased detection accuracy and leg tracking failure caused by the lidar sensor being placed on the ground.

Although sensors for measuring gait parameters have been proposed and evaluated in many studies, all literature on measuring gait parameters with lidar has adopted an offline model. In addition, in several algorithms based on other sensors, the fixed threshold (Li et al., 2021) and adaptive threshold methods (Bejarano et al., 2015) have been proposed.

In this paper, an online gait analysis algorithm applied to lower limb rehabilitation robot is designed. Based on the gait point cloud data collected by lidar, seven gait parameters including step length, stride length, gait cycle, stance time, step time, swing time and double support time are calculated with high accuracy in real time through four core algorithms: adaptive classification based on a sliding threshold, leg contour extraction, gait phase division and gait parameter calculation. Meanwhile, compared with the previous gait algorithm, which calculated the difference between the left and right legs in the y axis direction as the approximate spatial parameter, a spatial parameter-splicing algorithm based on a time series is proposed to effectively reduce the error influence of lidar coordinate changes and oblique walking trajectory on the gait spatiotemporal parameters. Based on the above gait algorithm, a dynamic window method is proposed, with 100 data points as the dynamic window. In addition, the method can identify the gait in real time at a sampling frequency of 10 Hz and calculate the gait parameters online. Finally, the optical motion capture system is used as a reference system to verify the feasibility and accuracy of the proposed gait analysis algorithm.

2 System overview

In this paper, a two-dimensional lidar sensor (RPLIDAR A2M8, Shanghai Silan Technology Co., LTD., Shanghai, China) with a scanning frequency of 5–15 Hz (typical value is 10 Hz) selected. The lidar sensor is connected with a driving board. External devices power and transmit data through two USB interfaces of the driving board, as shown in Fig. 1.

Figure 1 A photograph of the lidar sensor.

The robot used in this paper is the iReGo lower limb rehabilitation robot independently developed by our research group (J. C. Ji et al., 2018). The gait analysis system is divided into three parts: lidar, lower computer and upper computer of the lower limb rehabilitation robot. In order to more clearly describe the lower limb rehabilitation robot in this paper, Fig. 2 shows a 3D model drawing.

Figure 2 3D model of lower limb rehabilitation robot.

Lidar is used for the collection of leg data, and the data are transmitted to the upper computer of the lower limb rehabilitation robot through a USB cable. At the same time, the lidar sensor is powered by the terminal module of the lower computer system. The upper computer sends control instructions to the lidar sensor through the software interactive interface, so as to start and stop the sensor and stop the data transmission. In the meantime, the gait data from the lidar sensor are received for visual display and the real-time data and gait parameters are stored in the database for persistent management. The software interface is shown in Fig. 3.

Figure 3 The software interface of walking train.

Users need to wear a belt when using the robot for training. Figure 4 shows a scene of actually using the lower limb rehabilitation robot. Among them, the lidar sensor is installed in the control cabinet, which is sealed with a steel plate, and the control cabinet is slotted at the height of the laser emission light source so that the laser beam can be emitted from the control cabinet to obtain the user's gait information. Similarly, the relationship between the Cartesian coordinate system established by lidar and the human body can be clearly seen in Fig. 4. The opposite direction of human body walking is the positive direction of the y axis, and the measured angle of 0^∘ under polar coordinates is the positive direction of the x axis (i.e., the left direction of human body).

Figure 4 Practical scenario using the lower limb rehabilitation robot.

Since lidar can only detect a two-dimensional spatial plane, the installation height of the lidar sensor is adjusted to the height above the ground plane, and the distance between the tester's calf and the sensor is detected as the data basis of the gait parameter algorithm (see Fig. 5). Multiple tests in this paper have shown that the installation height should be 280 mm, which is designed to ensure that no shoes are detected during walking.

Figure 5 Influence of different lidar mounting heights on measurements.

3 Method

3.1 Algorithm

As shown in Fig. 6, the core of the gait analysis algorithm is to convert the original data collected by lidar into gait parameters, which is mainly composed of the following processes. Firstly, the point cloud data obtained from lidar scanning is preprocessed, and the range of data points is selected to eliminate the interference of miscellaneous points on the subsequent algorithm. Secondly, adaptive classification of selected data points, leg contour extraction and determination of left and right leg positions are carried out, in which the possible state of legs together is considered. In addition, according to the time series of the left and right leg positions in the walking process, the peak and trough points (namely toe off the ground, TO, and heel strike on the ground, HS) are identified. Finally, on the basis of gait phase division, the gait spatial parameter-splicing algorithm based on the time series is proposed to calculate the gait spatiotemporal parameters. The algorithm takes into account the variability of the lidar coordinate system applied to the robot, which can reduce the influence of its error on the gait spatiotemporal parameters. Based on the above gait analysis algorithm, 100 data points are used as the dynamic window to calculate the gait parameters of each frame and realize the real-time performance of gait analysis.

Figure 6 Flow chart of gait analysis algorithm.

3.2 Preprocessing

Define the arbitrary measurement data point p_ij obtained in the ith scan as

$$\begin{array}{}\text{(1)}& p_{ij}=\left(t_i,{\mathit{\theta }}_{ij},d_{ij},x_{ij},y_{ij}\right),\end{array}$$

where j is the index of the data points from each scan, t_i is the relative time at which the scan begins, θ_ij is the relative angle of the jth point in the ith scan, d_ij is the Euclidean distance between the jth point and the origin defined by lidar, and (x_ijy_ij) is the Cartesian coordinate converted from the original polar coordinate of the jth point.

Assume that the lidar scanning plane is parallel to the ground, its position with respect to the rehabilitation robot remains unchanged and the x axis is defined as the measured angle of 0 ^∘. The established Cartesian coordinate system is shown in Fig. 4, and the specific formula is

$$\begin{array}{}\text{(2)}& {\displaystyle }{x}_{ij}=d_{ij}\cos {\mathit{\theta }}_{ij},\text{(3)}& {\displaystyle }{y}_{ij}=d_{ij}\sin {\mathit{\theta }}_{ij}.\end{array}$$

Then the set of all laser points obtained in the ith scan is defined as follows:

$$\begin{array}{}\text{(4)}& S_i=\left\{p_{ij}\right\}j=\mathrm{1},\mathrm{2},\mathrm{\dots },n,\end{array}$$

where n is the number of points at the ith scan (namely 800 for RPLIDAR A2M8).

Since the position of the lidar sensor remains unchanged relative to the rehabilitation robot, in order to extract the data points corresponding to the leg contour, a rectangular measurement area is defined with a certain range of effectiveness. In this case, the set of laser points is defined as follows:

$$\begin{array}{}\text{(5)}& S_i^{\mathrm{W}}=\left\{p_{ij}\in S_i,X_{min}<x_{ij}<X_{max},Y_{min}<y_{ij}<Y_{max}\right\},\end{array}$$

where X_min, X_max, Y_min and Y_max are the limits of the measurement area in Cartesian coordinates, in this paper taken as $\left[-\mathrm{300},\mathrm{300},\mathrm{0},\mathrm{1200}\right]$.

3.3 Leg recognition

3.3.1 Adaptive classification based on sliding threshold

In order to detect leg position, clustering of discrete data points obtained by lidar is a key step in leg recognition. In this paper, an adaptive clustering algorithm based on a sliding threshold is adopted. When the distance between the current point  p_ij and the previous point $p_{i,j-\mathrm{1}}$ is less than a threshold value D_{ij_max}, it can be roughly classified as the same class of data points. Conversely, p_ij and $p_{i,j-\mathrm{1}}$ are in different categories: p_ij is the first point of the new category while $p_{i,j-\mathrm{1}}$ is the last point of the old category. It is specifically expressed as

$$\begin{array}{}\text{(6)}& D_{ij\mathit{\_}max}=d_{i,j-\mathrm{1}}\left({\displaystyle \frac{\sin \mathrm{\Delta }\mathit{\theta }}{\sin \left(\mathit{\lambda }-\mathrm{\Delta }\mathit{\theta }\right)}}\right)+\mathrm{3}\mathit{\sigma },\end{array}$$

where D_{ij_max} is the threshold calculated of the jth laser point in the ith scan, Δθ is the scanning angle resolution of the lidar sensor, σ is the measurement distance error and λ is the adjustment constant parameter. Since the angle resolution and measuring distance error of the lidar sensor itself are 0.7^∘ and 10 mm, after multiple experimental adjustments and verifications, in this paper, we take σ=0.01 m, λ=10^∘.

As shown in Fig. 7, the process of the proposed adaptive clustering method based on a sliding threshold can be clearly represented, where Δθ is amplified to simplify the visualization process. Figure 8 shows the classification results.

Figure 7 Schematic diagram of adaptive clustering principle.

Figure 8 Data classification of adaptive clustering based on a sliding threshold.

According to experimental experience, at least three points on each leg are required for leg circle fitting. If the number of data points in the category obtained after adaptive clustering is less than 3, such as some outliers, these data points will be excluded as environmental noise. There is also the possibility that the lidar sensor itself scans for less laser points.

Meanwhile, the distance between the jth laser point and the previous laser point in the ith scan is defined as

$$\begin{array}{}\text{(7)}& D_{ij}=\left|\left|p_{ij}-p_{i,j-\mathrm{1}}\right|\right|=\sqrt{{\left(x_{ij}-x_{i,j-\mathrm{1}}\right)}^{\mathrm{2}}+{\left(y_{ij}-y_{i,j-\mathrm{1}}\right)}^{\mathrm{2}}}.\end{array}$$

Then the clustering condition is

$$\begin{array}{}\text{(8)}& Q_i=\left\{\begin{array}{ll}{Q}_i+\mathrm{1}& D_{ij}>D_{ij\mathit{\_}max}\\ Q_i& D_{ij}\le D_{ij\mathit{\_}max}\end{array},\right.\end{array}$$

where Q_i is the number of categories of data points scanned at the ith time.

3.3.2 Leg contour extraction

The cross section between the laser scanning plane and the human leg can be approximated as a circle. Since the lidar sensor can only scan the contour of the calf from a certain angle range, the obtained point cloud data set is an approximately circular arc and the complete contour of the calf cannot be obtained, as shown in Fig. 8. Therefore, it is necessary to fit the circle as the approximate contour of the calf.

In order to take stability into consideration, this paper first uses a geometric method to calculate the matching circle, and determines the matching circle parameters according to the least square method, as follows:

$$\begin{array}{}\text{(9)}& min\sum _{j=\mathrm{1}}^n[{\left(x_{ij}-x\left({\mathit{\theta }}_{ij}\right)\right)}^{\mathrm{2}}+{\left(y_{ij}-y\left({\mathit{\theta }}_{ij}\right)\right)}^{\mathrm{2}}].\end{array}$$

Then the center and radius of the matching circle of laser data points in each frame can be obtained. However, through experimental tests, it is found that the real-time performance of this method cannot meet the predetermined requirements in the case of real-time data acquisition, so this paper simplifies it.

In the case of ideal clustering results (i.e., after adopting an adaptive clustering algorithm based on a sliding threshold, there are two classes of data points and the number of each class of data points reaches above the required number of center fitting), we assume a cluster m of ith scan, defining each data point in the cluster as $r_{ik}=\left\{x_{ik},y_{ik}\right\}$, $k=\mathrm{1},\mathrm{2},\mathrm{\dots },n_m$. Based on experimental experience, the circle with a fixed radius of R_leg was initially defined as the model for leg recognition. Estimate the location of the center point of each cluster based on geometric constraints:

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{u}_{im}& ={\displaystyle \frac{x_{i,ck}\times \left(\sqrt{x_{i,ck}^{\mathrm{2}}+y_{i,ck}^{\mathrm{2}}}+R_{\mathrm{leg}}\right)}{\sqrt{x_{i,ck}^{\mathrm{2}}+y_{i,ck}^{\mathrm{2}}}}}\\ & =\cos \mathit{\beta }\times \left(\sqrt{x_{i,ck}^{\mathrm{2}}+y_{i,ck}^{\mathrm{2}}}+R_{\mathrm{leg}}\right)\end{array},\text{(11)}& {\displaystyle }\begin{array}{rl}{w}_{im}& ={\displaystyle \frac{y_{i,ck}\times \left(\sqrt{x_{i,ck}^{\mathrm{2}}+y_{i,ck}^{\mathrm{2}}}+R_{\mathrm{leg}}\right)}{\sqrt{x_{i,ck}^{\mathrm{2}}+y_{i,ck}^{\mathrm{2}}}}}\\ & =\sin \mathit{\beta }\times \left(\sqrt{x_{i,ck}^{\mathrm{2}}+y_{i,ck}^{\mathrm{2}}}+R_{\mathrm{leg}}\right)\end{array},\end{array}$$

where ck is the intermediate value of the index of the data points contained in the various clusters.

$$\begin{array}{}\text{(12)}& ck=\text{median}\left\{k\right\},\end{array}$$

where R_leg is the fixed radius for fitting the center of the leg circle, which is taken as 40 mm in this paper considering the user's calf radius at the height of lidar scanning. $c_{im}=\left\{t_i,u_{im},w_{im}\right\}$ is the estimated position of the center point of the cluster in the ith scan. β is the angle between the origin and the central point of each cluster.

In this mathematical geometry model, the estimated center points of each cluster and the center points of the data point index are considered to be numerically approximate in the angle formed by the origin when the aspect ratio difference is large. Figure 9 shows a schematic diagram of leg center estimation in one scan, where the blue curve in the figure represents the right foot of the tester and the red curve represents the left foot of the tester. Figure 10 shows the actual fitting circle results.

Figure 9 Estimation of leg center position.

Figure 10 Contour fitting results.

3.3.3 Leg tracking and differentiation

When the clustering result is ideal, let $S_i^{\mathrm{WL}}$ represent the left leg data point cluster and $S_i^{\mathrm{WR}}$ represent the right leg data point cluster. On the premise that the center position of the two clusters is estimated, $L_{i,m,m}$ is defined as the Euclidean distance between the center point position of the category c_i,m obtained by the ith scan and the center point position of the previous scan $c_{i-\mathrm{1},m}$.

Thus, let the distance $L_{i,m,\mathrm{L}}$ of the left leg cluster and the distance $L_{i,m,\mathrm{R}}$ of the right leg cluster be defined as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{L}_{i,m,\mathrm{L}}=\left|\left|c_{i,m}-c_{i-\mathrm{1},\mathrm{L}}\right|\right|,\text{(14)}& {\displaystyle }{L}_{i,m,\mathrm{R}}=\left|\left|c_{i,m}-c_{i-\mathrm{1},\mathrm{R}}\right|\right|.\end{array}$$

Before left and right leg judgment, left and right leg cluster initialization is performed. According to the coordinate system and the fixed gait law, the position of the left and right legs corresponding to the two clusters are determined for the first time through the abscissa size of the center in the first frame, that is, the center point cluster with the smaller abscissa value is the right foot cluster.

From this, the following model can be obtained:

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{lll}{c}_{i,m}\in S_i^{\mathrm{WL}}& L_{i,m,\mathrm{L}}<L_{i,m{}^{\prime },\mathrm{L}}& \text{only left foot}\\ c_{i,m}\in S_i^{\mathrm{WR}}& L_{i,m,\mathrm{R}}<L_{i,m{}^{\prime },\mathrm{R}}& \text{only right foot}\\ c_{i,m}\in S_i^{\mathrm{WR}}& u_{im}<u_{im{}^{\prime }}& \text{both feet}\end{array}.\right.\end{array}$$

However, in the actual situation, clustering is often not ideal, as there may be more than two clusters that also meet the number of leg center fitting requirements. In this case, a judgment factor is proposed to select the first two categories with a larger value as data point clusters for leg fitting. The factor is defined as follows:

$$\begin{array}{}\text{(16)}& K={\displaystyle \frac{n_m}{L_{i,m,m}}},\end{array}$$

where n_m is the number of data points within cluster m obtained in the ith scan.

Considering the situation that legs may be together in the actual walking process, the threshold L is set when only one cluster exists after adaptive clustering. In this paper, L is 150 mm. When the distance between the first and the end points in the cluster is greater than the threshold value, it is considered that there are two types of clusters; otherwise, there is only one type of clusters. The results can be obtained according to the above model.

3.4 Gait feature detection based on peak values and geometric constraints

Human walking is a periodic process. It can be concluded that the leg coordinate time series of walking presents a curve similar to a sine function. In the actual process of walking, due to the influence of noise factors, the curve is prone to a tiny jitter phenomenon, which leads to pseudo-peaks and pseudo-troughs. This situation has great interference for the recognition of two gait key points: HS and TO. Based on this, a gait feature detection method based on peak values and geometric constraints is adopted in this paper. If the following two conditions are met, the feature point is considered effective.

•  

  Condition 1. Perform first-order curve fitting on the current gait data, excluding pseudo-troughs above the curve and pseudo-peaks below the curve.

•  

  Condition 2. For the same side, there is only one trough between adjacent peaks. For the opposite side, there is only one peak between adjacent peaks.

3.5 Calculation of gait parameters

After detecting the gait features of the point cloud data, it is represented by the smooth leg coordinate time series. The actual curve is shown in Fig. 11, which roughly presents a periodic law, but it is not as smooth as the ideal curve, with pseudo-troughs and pseudo-peaks.

Figure 11 Gait curve diagram.

According to the gait curve obtained, the peak value (TO) and trough value (HS) in the figure are extracted as eigenvalues. In this paper, HS time is defined as the time point at which the heel touches the ground, which is approximately the discrete time corresponding to the minimum value in the smooth time series. TO time is defined as the time point when the toe leaves the ground, which is approximately the discrete time corresponding to the maximum value in the smooth time series. HS and TO positions are defined as the distance from the lidar sensor in the forward and backward directions.

In this paper, seven spatiotemporal parameters are calculated (i.e., stance time, swing time, double support time, stride time, step time, stride length and step length). A diagram of gait time parameters is shown in Fig. 12. The formula is as follows:

Figure 12 Schematic diagram of gait time parameters.

In the kth gait cycle, T_HS(k) and T_TO(k) are assumed to be the corresponding moments of HS and TO, and $T_{\mathrm{TO}}^{\prime }\left(k\right)$ is assumed to be the moment that the opposite toe leaves the ground.

Taking the heel striking the ground as the starting point of a cycle, then gait cycle T is the time from one heel striking to the next heel striking on the same side, which can be expressed as

$$\begin{array}{}\text{(17)}& T=T_{\mathrm{HS}}(k+\mathrm{1})-T_{\mathrm{HS}}\left(k\right).\end{array}$$

Swing time T_sw is the time from the TO to the heel striking on the same side, expressed as

$$\begin{array}{}\text{(18)}& T_{\mathrm{sw}}=T_{\mathrm{HS}}\left(k\right)-T_{\mathrm{TO}}\left(k\right).\end{array}$$

Stance time T_st is the time from heel striking to the toe off on the same side, expressed as

$$\begin{array}{}\text{(19)}& T_{\mathrm{st}}=T_{\mathrm{TO}}\left(k\right)-T_{\mathrm{HS}}\left(k\right).\end{array}$$

Double support time T_ds is the time from heel striking to the opposite toe off, expressed as

$$\begin{array}{}\text{(20)}& T_{\mathrm{ds}}=T_{\mathrm{TO}}^{\prime }\left(k\right)-T_{\mathrm{HS}}\left(k\right).\end{array}$$

In the kth gait cycle, Y_HS(k) and $Y_{\mathrm{HS}}^{\prime }\left(k\right)$ are assumed to be the position coordinates when one heel strikes and the opposite heel strikes. The calculation method of standard gait spatial parameters is as follows:

$$\begin{array}{}\text{(21)}& \text{SL}=Y_{\mathrm{HS}}\left(k\right)-Y_{\mathrm{HS}}^{\prime }\left(k\right).\end{array}$$

Stride length is the longitudinal straight distance from one heel striking to the next heel striking on the same side, expressed as

$$\begin{array}{}\text{(22)}& \text{SR}=Y_{\mathrm{HS}}\left(k\right)-Y_{\mathrm{HS}}(k+\mathrm{1}).\end{array}$$

Since the measurement height of the lidar sensor is above the ankle height of the tester, and the shoe height worn by different people during the test also has an impact on the data measured by lidar. A height factor G is introduced based on the analysis of the human body model, as shown in Figs. 13 and 14.

Figure 13 Gait parameter relationship diagram based on human body model.

Figure 14 Gait parameter definition diagram.

According to the figure, the following two relationships can be obtained:

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{D_{\mathrm{1}}-D_{\mathrm{2}}}{L_{\mathrm{c}}}}={\displaystyle \frac{H_{\mathrm{P}}-H_{\mathrm{R}}}{H_{\mathrm{P}}}}\text{(24)}& {\displaystyle }{H}_{\mathrm{P}}^{\mathrm{2}}+{\left({\displaystyle \frac{L_{\mathrm{c}}}{\mathrm{2}}}\right)}^{\mathrm{2}}=L_{\mathrm{t}}^{\mathrm{2}},\end{array}$$

where L_c is the step length in a gait cycle and H_R and H_P are the height of the lidar sensor and the human hip joint from the ground, respectively. D₁ and D₂ are the distances from the tester's right and left legs to the origin, and L_t is the tester's leg length.

From this, the height factor can be obtained:

$$\begin{array}{}\text{(25)}& G={\displaystyle \frac{H_{\mathrm{P}}}{H_{\mathrm{P}}-H_{\mathrm{R}}}}\approx {\displaystyle \frac{L_{\mathrm{t}}}{L_{\mathrm{t}}-H_{\mathrm{R}}}}.\end{array}$$

Most previous research on lidar is based on the sensor being fixed to the ground plane to calculate the gait parameters. In this paper, the lidar sensor is applied to the rehabilitation robot. When calculating the gait parameters, especially the step length parameter with a different frame, the direction of the rehabilitation robot may not be straight forward during the walking process, and the walking trajectory of the tester may not be a straight line. Therefore, a spatial parameter-splicing algorithm based on a time series is proposed, as shown in Fig. 15.

Figure 15 Spatial parameter-splicing algorithm based on a time series.

Define the new center point of the leg position formed by splicing as $r_i^N$ in the ith scan, where $r_i^N=\left\{t_i,x_i^N,y_i^N\right\}$; then, the formulas for calculating the step length is as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle }{x}_i^N=x_i^{\mathrm{R}}-x_i^{\mathrm{L}}+x_{i-\mathrm{1}}^{\mathrm{L}},\text{(27)}& {\displaystyle }{y}_i^N=y_i^{\mathrm{R}}-y_i^{\mathrm{L}}+y_{i-\mathrm{1}}^{\mathrm{L}},\text{(28)}& {\displaystyle }k={\displaystyle \frac{y_i^N-y_{i-\mathrm{1}}^{\mathrm{R}}}{x_i^N-x_{i-\mathrm{1}}^{\mathrm{R}}}},\text{(29)}& {\displaystyle }{x}_i^T={\displaystyle \frac{\left(k^{\mathrm{2}}\cdot x_{i-\mathrm{1}}^{\mathrm{R}}+k\cdot \left(y_{i-\mathrm{1}}^{\mathrm{L}}-y_{i-\mathrm{1}}^{\mathrm{R}}\right)+x_{i-\mathrm{1}}^{\mathrm{L}}\right)}{k^{\mathrm{2}}+\mathrm{1}}},\text{(30)}& {\displaystyle }{y}_i^T=k\cdot \left(x_i^T-x_{i-\mathrm{1}}^{\mathrm{R}}\right)+y_{i-\mathrm{1}}^{\mathrm{R}},\end{array}$$

where a vertical line is drawn from the support foot to the track of the swing foot and the vertical point is defined as $t_i^T$, $t_i^T=\left\{t_i,x_i^T,y_i^T\right\}$. The Euclidean distance between $t_i^T$ and $r_i^N$ is defined as the step length.

Through this algorithm, the errors of gait space parameters measured in mobile robots (such as the following two cases) can be reduced, as shown in Figs. 16 and 17. In the difference algorithm (Sakdarakse and Somboon, 2020), the step length is the distance in the y direction between the legs at the HS point. Figure 16 shows a case where the walking track is oblique and the lidar coordinate system remains unchanged. sl and sw are the step length and step width calculated by the difference algorithm, and sl^′ and sw^′ are the step length and step width calculated by the algorithm proposed in this paper. Similarly, the stride length is the distance on the diagonal line.

Figure 16 Walking path oblique condition.

Figure 17 Lidar coordinate system changes.

Figure 17 shows a situation when the track is straight and the lidar coordinate system changes. For the angle change of the lidar sensor θ, a=a₁, sl, sl₂ and sl₃ are the actual step length, the step length calculated by the proposed algorithm and the step length calculated by the difference algorithm, respectively. As can be seen from the figure, the parameters calculated by the proposed algorithm in this paper are closer to the actual parameters, and as θ increases and the spatial parameters in actual walking are larger, the error can be effectively reduced by this algorithm.

3.6 Online real-time analysis algorithm

Based on the above gait analysis algorithm, the post-calculation of each gait parameter can be realized in the case of data sets. However, the real-time ability of the algorithm to recognize heel striking should be good. In this paper, 100 data points are taken as a dynamic window and 10 Hz is taken as a sampling frequency. In addition, on the basis of the 100 data points, whether the 99th data point is a trough (that is, whether the leg heel strikes the ground) is identified. Through this dynamic window, the function of real-time recognition of gait and calculation of gait parameters is realized.

Finally, through the communication between the upper and lower computers, the gait of the patient is collected in real time and presented on the interface of the rehabilitation robot, which brings the fun of walking training to the patient. It also provides gait parameter data for doctors to evaluate patients' recovery status, and the data are synchronously saved to the cloud database so that doctors can remotely view patient reports at any time.

4 Experiment and results

4.1 Experimental setup

The optical motion capture system (OptiTrack, NaturalPoint, Inc, Corvallis, USA) is used as the reference system. In order to determine the joint location, optical markers are placed on the testers' knee joints and lateral ankles. A total of 22 markers are attached to each tester's lower limbs (see Fig. 18 for detailed marker locations). HS time (heel strike time) and TO time (toe off time) are automatically estimated based on the marked motion path. HS time is defined as the time when the distance between the marks of the left and right lateral malleolus of the leg maximizes in the direction of the tester (Tanabe et al., 2017). The TO time is defined as the maximum angle between the calf (the line connecting the center mark of the knee joint rigid body and the center mark of the ankle joint rigid body) and the vertical line in the sagittal plane. Based on this, HS position and TO position are calculated according to time. Therein, the origin of the coordinate system set by the reference system is located at the middle point of the walking distance.

Figure 18 Wearing details of the experimental tester.

The OptiTrack Prime 17W camera is set around the top of the laboratory (see Fig. 19, where individual motion capture cameras are represented as M) to form the capture area of the entire motion capture system. The rehabilitation robot equipped with lidar is located at one end of the test area (represented by T in Fig. 19).

Figure 19 Experimental environment settings.

4.2 Verifying the effectiveness of the gait analysis algorithm

To verify the validity of the proposed gait analysis algorithm, six young testers (three males and three females: average age 24–26, height 170–175 cm and weight 50–60 kg) are recruited to participate in this experiment. In this experiment, each tester carries out five experiments. In the scene of using the rehabilitation robot, gait parameters calculated by the proposed gait analysis algorithm are compared with those measured by the OptiTrack.

The experiment scenario is as follows: before the experiment, the leg length and height of the shoes worn by the tester are recorded, and the tester wears the belt of the rehabilitation robot. The robot moves forward with the force sensor in the belt, and there is a certain distance between the tester's waist and the robot. In the five experiments, the first three experiments are carried out according to the fixed step length and the number of steps marked on the ground (the fixed step length is 350 mm according to the height of the tester in this experiment). For the last two times, it is enough to follow a certain length of a straight path according to the usual walking habits of the tester, taking into account the differences in walking speed and stride length among the different individuals. At the same time, the optical markers are affixed to 22 positions of the lower limbs of the tester. During the experiment, the synchronization of the gait analysis algorithm and data collected by 3D motion capture system is ensured.

4.3 Analysis of gait parameter results

The statistics and errors of gait parameters detected by the reference system and the proposed gait analysis algorithm are shown in Table 1.

Table 1 Comparison of spatiotemporal parameters between the reference system and the proposed gait analysis algorithm.

As can be seen from Table 1, for the measurement of gait parameters, because the values of swing time and double support time are relatively small compared with other time parameters, the relative errors of the two are relatively large, 23.87 % and 14.27 %, respectively. The relative errors of the other parameters are all kept below 8 %. At the same time, from the measured step length and step width, the measured values of 360.3, 714.2, 368.3 and 728.6 mm in the experimental group and the control group are slightly different from the fixed step length of 350 mm and step width of 700 mm specified in the experiment. In summary, compared with other gait analysis systems (Fudickar et al., 2018; Iwai et al., 2019), the measurement error of lidar gait analysis algorithm is within the acceptable range. In addition, it can be observed that the walking cycle time of healthy adults is about 1.3 s, the swing time is about 0.47 s and the stance time is about 0.83 s (Kotiadis et al., 2010). On the lower limb rehabilitation robot, the corresponding time parameter value is significantly extended, indicating that rehabilitation walking in the lower limb rehabilitation robot is different from walking in natural, unrestrained conditions.

The correlation and consistency between the two systems can be evaluated by the Pearson correlation coefficient and intra-group correlation coefficient (ICC(2,1)) and 95 % confidence interval. When using Pearson correlation analysis, it is necessary to judge its significance (p) and correlation coefficient (r). Significance indicates whether there is a relationship between the two system, while the correlation coefficient is used to judge the strength of the relationship between the two systems. The correlation of gait parameters between the two systems is shown in Fig. 20. Table 2 shows the correlation and consistency of gait parameters measured between the two systems.

Figure 20 Correlation of gait parameters between the two systems: (a) gait cycle, (b) double support time, (c) stance time, (d) swing time, (e) step time, (f) step length and (g) stride length.

Table 2 Correlation of gait parameters between the proposed system and the reference system.

It can be seen from Fig. 20 and Table 2 that all parameters present a significant positive correlation (p<0.01), among which all gait parameters except swing time show good correlation and consistency with the dynamic arrest system (r>0.75, ICC(2,1) >0.75). The gait cycle and the stance time show a very high correlation and consistency (r=0.99, ICC(2,1) =0.99 and r=0.98, ICC(2,1) =0.98), while the swing time has relatively low correlation and consistency compared with other time parameters (r=0.47, ICC(2,1) =0.45). This may be caused by two factors. First, the swing time parameter itself has a small measurement range, and ICC is the ratio of the true variance to the true variance plus error. If the true variance is small, the ICC will be relatively small. Secondly, there are certain errors in the measurement method itself. The proposed algorithm will have similar errors no matter how the measurement parameters change. Therefore, the decrease of the measured value with the measurement error may reduce the concurrent validity of the time parameter measured in the proposed method (Iwai et al., 2019).

Bland–Altman (BA) diagram is a common clinical method used to evaluate the consistency of two measurement methods. At the same time, in order to prove the influence of individuals in each gait parameter on the difference, a BA diagram of each gait parameter is drawn for this paper (as shown in Fig. 21). The horizontal axis represents the mean value of the gait parameters of the two systems. The vertical axis represents the difference between the measured values of the two systems. The blue horizontal solid line in the middle represents the mean value of the difference. The orange horizontal dashed line represents the position where the mean value of the difference is 0. The upper and lower brown horizontal dashed lines respectively represent the upper and lower limits of 95 % consistency limits, namely, the standard deviation of 1.96 times.

Figure 21 Bland–Altman diagrams of the two systems with respect to gait parameters: (a) gait cycle, (b) double support time, (c) stance time, (d) swing time, (e) step time, (f) step length and (g) stride length.

If the difference of the measurement results follows a normal distribution, then 95 % of the difference should lie within mean ± 1.96 SD. When most of the difference points lie within this interval and the mean value of difference is close to zero, this indicates that the consistency of the two systems is better. As can be seen from Fig. 21, most of the points fall within the consistency limit and the average value of the difference is close to zero, which proves that the proposed gait analysis algorithm and the reference system are very consistent.

4.4 Result analysis of different spatial parameter calculation methods

For the same point cloud data, the proposed spatial parameter-splicing algorithm based on the time series and the existing difference algorithm are used to compare step length results. The statistical values and errors obtained are shown in Table 3.

Table 3 Comparison of step length of the two algorithms.

As can be seen from Table 3, when the tester walks randomly along a linear path of a certain length in accordance with their usual walking habits, the relative error and absolute error of the step length calculated by the proposed algorithm are 4 % and 14.2 mm. Compared with the relative error of 5.9 % and the absolute error of 21.2 mm calculated by the existing difference algorithm, it can be seen that the error calculated by the proposed algorithm is smaller than that calculated by the difference spatial parameter algorithm.

4.5 Analysis of the real-time effect of the proposed algorithm

In the above experiment, when the tester used the lower limb rehabilitation robot for walking training, the robot interactive interface simultaneously showed the footprints of the tester. For example, the tester steps out on their left leg and we record the moment when their left leg reaches HS. The proposed gait analysis algorithm synchronously detects the data of the left leg. The processed data are transmitted to the robot through the interface based on TCP communication, and the footprint image and parameters are synchronously displayed on the interface.

Table 4 Average recognition delay and RMSE of the proposed algorithm.

As can be seen from Table 4, the average recognition delay of the proposed gait analysis algorithm is about 0.486 s, which is within the acceptable range. The recognition delay is mainly the detection delay of an HS event and the fluctuation of the leg data in the HS event. The scanning frequency of the lidar sensor itself is 10 Hz. According to the law of the gait curve, after reaching the HS point, the HS situation can be determined only when the next data point is detected within at least 0.1 s.

5 Conclusion

At present, most quantitative gait analysis is usually limited by site and scope and cannot be well integrated with lower limb rehabilitation robots. In this paper, a gait analysis algorithm for a lower limb rehabilitation robot is proposed, which uses lidar to obtain gait phase information and gait spatiotemporal parameters. The algorithm includes the adaptive clustering method based on a sliding threshold, leg contour fitting and center estimation, recognition of gait phase key points and calculation of gait time and space parameters. To effectively reduce the errors caused by patient trajectory deviation and changes in the lidar coordinate system due to robot motion during walking, a spatiotemporal parameter-splicing algorithm based on a time series is proposed. Compared with the gait parameters measured by the optical motion capture system, the experimental results verify the effectiveness and accuracy of the proposed algorithm. Except for the two parameters of double support time and swing time, which have relatively large relative errors due to their small values, the relative errors of other gait parameters are kept below 8 %. Based on the gait analysis algorithm, a dynamic window method is proposed in this study, which takes 10 Hz as the sampling frequency and 100 data points as the dynamic window to identify whether the 99th data point is the HS point, so as to realize real-time recognition and gait parameter calculation. The experimental results show that the average recognition delay of the proposed algorithm is about 0.486 s, which is within the acceptable range, and it shows good performance for the application of online gait parameter calculation of mobile robots.

In future development work, emphasis should be placed on the application of online gait parameters in real-time robot control and interaction strategies, intelligent recognition of hemiplegia and other abnormal gaits, adaptive adjustment of rehabilitation plan, etc., which have good application prospects in the fields of mobile robots, rehabilitation and diagnostic monitoring.