A complete process for human dynamics analysis including parameter identification and sEMG-torque estimation

Abstract

In order to analyze the motion of the human lower limb, a comprehensive process was created in this paper. Firstly, an exoskeleton robot platform was leveraged and a dynamic model of human-exoskeleton lower limb was created for identifying human parameters. Then the sEMG signal was utilized to extract information about muscle activity through multi-group experiments, and a backpropagation neural network (BPNN) was built to forecast joint torque. An inverse dynamics analysis combining the human motion data with the dynamic model can not only verify the reliability of prediction result by this BPNN but also the correctness of identified results before. Moreover, the mean absolute error (MAE), root mean square error (RMSE), and Pearson correlation coefficient (PCC) were used as evaluation index of both identification and prediction results. The proposed protocol can give accurate identified parameters for subject and estimated joint torque from sEMG during swing motion. We believe it can be extended to various types of human motion movement and potentially applied to complete human motion analysis.

Introduction

Exoskeletons for the lower limbs are employed in the civilian sector to lighten the load on the human body and in the medical rehabilitation sector as assistive rehabilitation tools to help patients recover. ¹ Lower-limb exoskeletons aim to replicate human movements using mechanical components,^2,3 but their performance is constrained by the fact that joint moments of human are continuously and unconsciously modulated while moving around. The study of human motion is essential, but a significant obstacle is predicting joint moments in order to fully use the huge potential of lower limb exoskeletons. ⁴ The optimization and development of lower limb exoskeletons depend heavily on joint torque estimate, sometimes referred to as kinematic intent. Real-time estimation of muscle forces in the human lower extremity is a valuable assessment metric. ⁵ By observing the muscle forces produced by patients during rehabilitation, therapists can assess whether their treatment is within safe bounds, and patients can assess whether they are producing enough force to improve their physical condition. ⁶ However, it is currently difficult to estimate joint moments in the lower extremities in real time during human locomotion. ⁷ This is because it is challenging to estimate lower extremity moments due to the existing actuator redundancy, many nonlinearities, intrusive experimental settings, and ethical issues.

It is a difficult task to choose the right signals to forecast the user’s motor intention. ⁸ The majority of current clinical work and research focuses on the prediction of moments using methods including direct analysis of electromyography (EMG) patterns, ⁹ extracted EMG envelopes, ¹⁰ and EMG driver models, ¹¹ together with joint moment analysis using dynamometer or inverse kinematic data. ¹² The surface electromyography (sEMG) signal is a thorough superposition of neural control signals from human surface muscles that can represent the state of contraction in human muscles. ¹³ The signal generation time is around 20–80 ms earlier than the time at which the muscle contracts. ¹⁴ Therefore, sEMG signals can be used to fully estimate the movement intention without causing any information loss or delay. The user’s purpose for various walking motions may be clearly seen in sEMG signals since they are non-intrusive, instructive, and pose no risk to the user. ¹⁵ Moreover, sEMG signals exhibit shorter time delay and a greater signal-to-noise ratio in current detection techniques than mechanical torque signals. ¹⁶ In conclusion, surface EMG signals are generally considered to be the best source for foretelling active joint moments. ¹⁷ In order to simulate the pattern identification of the acquired sEMG signals with joint torques, researchers now employ machine learning algorithms, such as neural networks and regression approaches. ¹⁸ By extracting sEMG signals from several distinct trials, this enables the prediction of joint moments by identifying complicated links between the inputs (sEMG signals) and outputs (joint torques). ¹⁹

Currently, sensor acquisition can be used to get sEMG signals, but joint torques can only be acquired indirectly. ²⁰ Inverse dynamics analysis is necessary to estimate the torques produced by the lower limb muscle groups from a mechanical perspective, ²¹ which is called the gold standard in human gait analysis. ²² The necessary inertial parameters, such as the mass, center of mass location, and moment of inertia of the human lower limb, must be accurately estimated in order to calculate the torques produced by the lower limb muscle groups using inverse dynamics. ²³ It is challenging to determine the inertial parameters of lower limb exoskeletons since they are mostly made of mechanically structured electronic modules that are controlled by motors or electro-hydraulics with various control systems. ²⁴ Due to the variations in each individual user’s physical characteristics, the calibration of the inertial parameters to the particular user is required in order to be able to forecast outcomes for various users. ²⁵ The accuracy of the inertial parameters of the human lower limbs is crucial for determining the joint torques obtained through inverse dynamics calculations, so it is necessary to identify and separate the inertial parameters of the exoskeleton and the human lower limbs.

In this paper, dynamic models of human lower limbs and human-exoskeleton coupling were established at first. The human inertia parameters were then separated from the human-exoskeleton coupling model through kinetic analysis and rational experiments. Multiple sets of sEMG-joint torque data were gathered in order to create the neural network’s subset. In this instance, the mapping of sEMG data to joint torques was predicted using a backpropagation neural network (BPNN). The dynamic model was paired with motion data, such as the joint angle, to carry out the suggested dynamics analysis, which may confirm the correctness of the BPNN prediction and the parameter identification.

Methods

After being fully informed about the experimental procedure, a 1.7 m/50 kg healthy female subject signed her informed consent to participate in the study, which was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the University of Electronic Science and Technology of China (106142023021725090). The complete process for human dynamics analysis proposed in this paper consists of three main parts as Figure 1: (1) Experiment (a) for dynamic parameters identification for human lower limbs; (2) Experiment (b) for Neural network training for predicting human joint torque from sEMG signals; and (3) Experiment (c) for the trained neural network validation.

Figure 1.

(a–c) The three experiments of process.

Dynamic parameters identification for human lower limbs

Figure 2 shows the schematic diagram of a 2-DOF lower extremity exoskeleton coupled with human lower limb, where the corresponding symbols are listed in Table 1. As the Lagrangian function, which is the difference between the kinetic and potential energy of a mechanical system, takes simply computational effort, ²⁶ the dynamic model of human lower limbs is constructed using Lagrange modeling approach as follows ²⁷ :

$${\mathrm{M}}_{\mathrm{h}}\left(\mathit{\theta }\right)\stackrel{··}{\mathit{\theta }}+{\mathrm{C}}_{\mathrm{h}}(\mathit{\theta },\stackrel{·}{\mathit{\theta }})\stackrel{·}{\mathit{\theta }}+{\mathrm{G}}_{\mathrm{h}}\left(\mathit{\theta }\right)={\mathbf{\tau }}_{\mathrm{r}}+{\mathbf{\tau }}_{\mathrm{h}}+{\mathbf{\tau }}_{\mathrm{ext}},$$

(1)

where $\mathit{\theta }=[{\theta }_1,{\theta }_2{]}^T\in {\mathrm{R}}^2$ represents the joint angle vector of human and exoskeleton’s lower limbs, the dot over $\mathit{\theta }$ represents the derivative with respect to time, ${\mathrm{M}}_{\mathrm{h}}\left(\mathit{\theta }\right)$ , ${\mathrm{C}}_{\mathrm{h}}(\mathit{\theta },\stackrel{·}{\mathit{\theta }})$ , and ${\mathrm{G}}_{\mathrm{h}}\left(\mathit{\theta }\right)$ represent the Inertia, Coriolis, and Gravity matrices of human lower limbs, ${\mathbf{\tau }}_{\mathrm{h}}\in {\mathrm{R}}^2$ and ${\mathbf{\tau }}_{\mathrm{ext}}\in {\mathrm{R}}^2$ represent joint torques from human muscle contraction and external environment (the exoskeleton). In a similar manner, we also obtain the governing equation of the exoskeleton as follows:

$${\mathrm{M}}_{\mathrm{r}}\left(\mathit{\theta }\right)\stackrel{··}{\mathit{\theta }}+{\mathrm{C}}_{\mathrm{r}}(\mathit{\theta }\mathbf{,}\stackrel{·}{\mathit{\theta }})\stackrel{··}{\mathit{\theta }}+{\mathrm{G}}_{\mathrm{r}}\left(\mathit{\theta }\right)={\mathbf{\tau }}_{\mathrm{r}}+{\mathbf{\tau }}_{\mathrm{f}}-{\mathbf{\tau }}_{\mathrm{ext}},$$

(2)

where ${\mathbf{\tau }}_{\mathrm{r}}\in {\mathrm{R}}^2$ denotes the vector of exoskeleton actuator torque applied at the joints, ${\mathbf{\tau }}_{\mathrm{f}}\in {\mathrm{R}}^2$ denotes the exoskeleton friction torque vector, ${\mathrm{M}}_{\mathrm{r}}\left(\mathit{\theta }\right)$ , ${\mathrm{C}}_{\mathrm{r}}(\mathit{\theta }\mathbf{,}\stackrel{·}{\mathit{\theta }})$ , and ${\mathrm{C}}_{\mathrm{r}}\left(\mathit{\theta }\right)$ represents the Inertia, Coriolis, and Gravity matrices of exoskeleton robot respectively. Adding the two equations together yields the dynamic equation of the coupled system as follows ²⁸ :

$${\mathrm{M}}_{\mathrm{hr}}\left(\mathit{\theta }\right)\stackrel{··}{\mathit{\theta }}+{\mathrm{C}}_{\mathrm{hr}}(\mathit{\theta },\stackrel{·}{\mathit{\theta }})\stackrel{··}{\mathit{\theta }}+{\mathrm{G}}_{\mathrm{hr}}\left(\mathit{\theta }\right)={\mathbf{\tau }}_{\mathrm{hr}}+{\mathbf{\tau }}_{\mathrm{f}},$$

(3)

where ${\mathrm{M}}_{\mathrm{hr}}={\mathrm{M}}_{\mathrm{h}}+{\mathrm{M}}_{\mathrm{r}}$ , ${\mathrm{C}}_{\mathrm{hr}}={\mathrm{C}}_{\mathrm{h}}+{\mathrm{C}}_{\mathrm{r}}$ , ${\mathrm{G}}_{\mathrm{hr}}={\mathrm{G}}_{\mathrm{h}}+{\mathrm{G}}_{\mathrm{r}}$ , ${\mathbf{\tau }}_{\mathrm{hr}}={\mathbf{\tau }}_{\mathrm{h}}+{\mathbf{\tau }}_{\mathrm{r}}$ . Notably, we only construct the dynamic model in this study as the attention is only paid to motion experiments, the kinematic model such as Ayiz and Kucuk ²⁹ which is meaningful for the control system is not considered.

Figure 2.

Coupled model of lower extremity exoskeleton and human body.

Table 1.

Symbolic representation of physical parameters in Figure 2.

Symbol	Quantity	Unit
${\theta }_1/{\theta }_2$	Hip/knee: joint angle	rad
$L_1/L_2$	Thigh/shank: length	m
$c_1/c_2$	Thigh/shank: center of mass	m

As no sensors can be used to measure human joint torques, ${\mathbf{\tau }}_{\mathrm{h}}$ , it is impossible to identify ${\mathrm{M}}_{\mathrm{h}}$ , ${\mathrm{C}}_{\mathrm{h}}$ , and ${\mathrm{G}}_{\mathrm{h}}$ directly by using equation (1). Therefore, we propose a three-step procedure for the identification of dynamic parameters of human lower limbs. The first step uses equation (2) and a swing motion of the exoskeleton without a wearer $({\mathbf{\tau }}_{\mathrm{ext}}=0)$ to identify the parameters of the exoskeleton, that is, ${\mathrm{M}}_{\mathrm{r}}$ , ${\mathrm{C}}_{\mathrm{r}}$ , ${\mathrm{G}}_{\mathrm{r}}$ , and ${\tau }_{\mathrm{f}}$ . The second step uses equation (3) and a coupling dynamic experiment without human muscle contraction $({\mathbf{\tau }}_{\mathrm{h}}=0)$ to identify ${\mathrm{M}}_{\mathrm{hr}}$ , $C_{\mathrm{hr}}$ , and ${\mathrm{G}}_{\mathrm{hr}}$ . Finally, human parameters are obtained by using ${\mathrm{M}}_{\mathrm{h}}={\mathrm{M}}_{\mathrm{hr}}-{\mathrm{M}}_{\mathrm{r}}$ , ${\mathrm{C}}_{\mathrm{h}}={\mathrm{C}}_{\mathrm{hr}}-{\mathrm{C}}_{\mathrm{r}}$ , ${\mathrm{G}}_{\mathrm{h}}={\mathrm{G}}_{\mathrm{hr}}-{\mathrm{G}}_{\mathrm{r}}$ in the last step.

In the first and second steps, least square method is adopted for the identification with measure robot joint torques and angles, ${\mathbf{\tau }}_f$ and $\mathit{\theta }$ . To this end, equations (2) and (3) are converted into the following linear form:

$$A_{\mathrm{i}}(\mathit{\theta },\stackrel{·}{\mathit{\theta }},\stackrel{··}{\mathit{\theta }})X_{\mathrm{i}}={\mathbf{\tau }}_{\mathrm{i}},$$

(4)

where $\mathrm{i}$ is the index of the first $(\mathrm{i}=\mathrm{r})$ and second $(\mathrm{i}=\mathrm{hr})$ steps, $X_{\mathrm{i}}\in {\mathrm{R}}^8$ represents the parameter vector to be identified, ${\mathrm{A}}_{\mathrm{i}}(\mathit{\theta },\stackrel{·}{\mathit{\theta }},\stackrel{··}{\mathit{\theta }})\in {\mathrm{R}}^{2\times 8}$ represents regression matrix composed of joint angles, angular velocities and accelerations, and ${\mathbf{\tau }}_{\mathrm{i}}\in {\mathrm{R}}^{2\times 1}$ represents the joint torque. Details can be found in Appendix.

Excitation trajectory design: To improve the identification accuracy, sampling noise should be suppressed by designing a suitable excitation trajectory. Though the cubic splines ³⁰ and polynomial splines ³¹ have been adopted in many trajectory planning researches, considering the periodicity of motion, the $n$ -th order Fourier series is used to realize a general representation of the excitation trajectory as follows ³² :

$${\mathit{\theta }}_{\mathrm{e}}=\sum _{k=0}^n(a_k\sin \left(k\omega t\right)+b_k\cos \left(k\omega t\right)),$$

(5)

where $\omega$ is the excitation frequency, ${\mathrm{a}}_k$ and ${\mathrm{b}}_k$ are the coefficient vectors to be optimized. The algorithm used to getting optimal ${\mathrm{a}}_k$ and ${\mathrm{b}}_k$ in this paper is particle swarm optimization (PSO), which designs particles only with speed and positional properties. Each particle applies random initialization in the prescribed search space, with capabilities retaining the optimal position. Based on the swarm’s best location, the properties of particles are modified, which is similar to simulating birds in a flock.^33,34

Since the excitation trajectory should be related to the identification effect, the objective function of PSO in this paper was taken as the condition number of the regression matrix $A_{\mathrm{i}}$ . Besides, to ensure the safety both of the exoskeleton system and human, the constraint intervals for the excitation trajectory should be added, which considering the design threshold of the exoskeleton device and the joint angle limitation range during human movement ³⁵ are listed in Table 2. In summary, the optimization problem can be reduced to the following standard form:

$$\begin{array}{c}\min \mathrm{Fit}({\mathrm{A}}_{\mathrm{i}})=\mathrm{Cond}({\mathrm{A}}_{\mathrm{i}})\hfill \\ \mathrm{s}.\mathrm{t}.\mathit{\theta },\stackrel{·}{\mathit{\theta }},\stackrel{··}{\mathit{\theta }}\in \mathrm{Table}2\hfill \end{array}$$

(6)

Table 2.

The constraint intervals for the excitation trajectory.

Constraint terms	Constraint intervals	Unit
Hip joint angle	[0.0872, 1.2217]	rad
Knee joint angle	[−1.6232, −0.2269]	rad
Hip/knee angular velocity	[−1.6449, 1.6449]	rad/s
Hip/knee angular acceleration	[−5.1677, 5.1677]	rad/s2

Specifically speaking, the coefficient vectors and ${\mathrm{b}}_k$ were set as random particles with position and velocity range of [−50, 50] and [−3, 3] respectively in the PSO algorithm, where the particle swarm parameters containing swarm size of 60, learning factors of 1.1 and 1.7, inertial weights within [0.1, 0.9] and the maximum number of iterations of 1000 were considered. The objective function $\mathrm{Fit}({\mathrm{A}}_{\mathrm{i}})$ would be minimized until the iterations up to 1000.

Experimental data collection: For human lower limbs parameters identification, the data was collected at the exoskeleton robot as Figure 1(a) with the frequency of 100 Hz. ³⁶ The primary construction of the exoskeleton robot was built of steel, and two servo motors (GDM1-100N2/120N2) were used to move the hip and knee joints, which could be monitored by two absolute encoders (INC-4-150/INC-3-125), while the forces that interact with the human were measured by four 3D force sensors (JNSH-2-10kg-BSQ-12). Throughout the whole experiment, the subjects did not feel any discomfort. At the first step of the experiments for exoskeleton robot identification, the exoskeleton swung by itself after entering the desired excitation trajectory, which was generated by PSO mentioned before. At the second step for human body identification, the subject was standing connected with the exoskeleton during the trial, enabling the entire leg to swing and the knee to be in a natural posture. ³⁷ The swing data set was then established by gathering the real-time angle, angular velocity, angular acceleration, and joint torque information from exoskeleton system.

Assuming that there are $N$ sets of sampling points in the experiment, the sampling regression matrix ${\overline{\mathrm{A}}}_{\mathrm{i}}\in {\mathrm{R}}^{2N\times 8}$ and the exoskeleton torque vector ${\overline{\mathbf{\tau }}}_{\mathrm{r}1}/{\overline{\mathbf{\tau }}}_{\mathrm{r}2}\in {\mathrm{R}}^{2N\times 1}$ in the first and second experiment step can be obtained. The parameter vectors of exoskeleton to be identified in the fist step can be expressed as:

$${\mathrm{X}}_{\mathrm{r}}=({{\overline{\mathrm{A}}}_{\mathrm{r}}}^{\mathrm{T}}{\overline{\mathrm{A}}}_{\mathrm{r}}{)}^{-1}{{\overline{\mathrm{A}}}_{\mathrm{r}}}^{\mathrm{T}}{\overline{\mathbf{\tau }}}_{\mathrm{r}1},$$

(7)

Then the identified friction torque vector ${\hat{\tau }}_{\mathrm{f}}$ can be obtained using $X_{\mathrm{r}}$ , the parameter vectors of human-exo to be identified in the second step can therefore be obtained:

$${\mathrm{X}}_{\mathrm{hr}}={\left({{\overline{\mathrm{A}}}^{\mathrm{T}}}_{\mathrm{hr}}{\overline{\mathrm{A}}}_{\mathrm{hr}}\right)}^{-1}{{\overline{\mathrm{A}}}^{\mathrm{T}}}_{\mathrm{hr}}({\overline{\mathbf{\tau }}}_{\mathrm{r}2}-{\overline{\mathbf{\tau }}}_{\mathrm{f}}).$$

(8)

The identification quality can be attained by comparing the actual joint torque ${\overline{\mathbf{\tau }}}_{\mathrm{r}1}/{\overline{\mathbf{\tau }}}_{\mathrm{r}2}$ and the computed joint torque ${\hat{\mathbf{\tau }}}_{\mathrm{r}1}/{\hat{\mathbf{\tau }}}_{\mathrm{r}2}$ using identified parameters, which can be written as:

$$\begin{array}{c}{\hat{\mathbf{\tau }}}_{\mathrm{r}1}={\overline{\mathrm{A}}}_{\mathrm{r}}{\hat{\mathrm{X}}}_{\mathrm{r}},\hfill \\ {\hat{\mathbf{\tau }}}_{\mathrm{r}2}={\overline{\mathrm{A}}}_{\mathrm{hr}}{\hat{\mathrm{X}}}_{\mathrm{r}}-{\hat{\mathbf{\tau }}}_{\mathrm{f}}.\hfill \end{array}$$

(9)

Neural network training for predicting human joint torque from sEMG signals

Muscle activation is a crucial indicator of human motion, as well as the dynamic information including torque. Surface electromyography (sEMG) can investigate muscle function through the inquiry of myoelectric signals, which is the electrical signal the muscles emanate. ³⁸ To obtain more precise information on human motion, it’s necessary to establish sEMG-joint torque model. A backpropagation neural network (BPNN) was here trained to predict the lower limbs joint torque using sEMG in order to avoid complicated musculoskeletal modeling computations.

sEMG signal feature extraction: For sEMG-Torque data acquisition, each muscle has the capacity to influence several joints, and each joint requires the activity of multiple agonist and antagonist muscle groups. ³⁹ In this paper, four muscles shown in Figure 3– the rectus femoris (RF), vastus lateralis (VL), vastus medialis (VM), and biceps femoris (BF) – were chosen for sEMG collection with considering relevant joint torques in swing motion. Before experimental trials, maximum voluntary contraction (MVC) exercises were conducted for each of the four muscles.

Figure 3.

Muscles used to acquire sEMG signals.

For standardization, the sEMG signals were distributed among the maximal voluntary contraction (MVC) of appropriate muscles. The Fast Fourier Transform (FFT) was utilized to extract five time features of sEMG signals, including mean value, root mean square (RMS), peak value, variance, and numbers of zero crossing points as well as three frequency features including gravity frequent, RMS frequency, and frequency standard deviation. All data were exported to MATLAB for analysis.

BPNN construction: For the human lower limbs swing model constructed in this paper, the torques at the hip and knee joints are estimated with corresponding sEMG signals using two BP neural networks respectively. The simplest BPNN structure of only three layers was conducted based on Neural Network Toolbox of MATLAB, with “tansig” and “purelin” as the activation function for hidden and output layers. The input data matrix is all extracted features of sEMG signals mentioned above, which has 32 nodes for input layer and the output is the joint torque with only 1 node in this paper. As for the number of hidden layer nodes, we first calculated the range of its values by ⁴⁰ :

$$n_h=\sqrt{n_i+n_o}+a,$$

(10)

where $n_{\mathrm{i}}$ is the number of nodes in the input layer, $n_{\mathrm{o}}$ denotes the number of nodes in the output layer, and $a$ is an integer between 30 and 40 in this paper. The exact number of hidden layer was decided according to the lowest mean square error (MSE) of the network performance during each value testing in training.

As long as the structure of BPNN is determined, the network can be trained by inputting training data subset from the above experiment. The raw data was automatically divided into three parts as training set for 70%, validation set for 15%, and test set for 15% using “Random” algorithm in the toolbox. The Levenberg-Marquardt optimization ⁴¹ was used with initial damping factor $\mu =0.001$ . The number of training iterations, learning rate, and minimum error of target was set as 1000, 0.01, and $1\times {10}^{-6}$ , respectively.

Experimental data collection: The data for network training were acquired at the stiff platform as Figure 1(b), where a 3D force sensor (LH-SZ-02-100N) was used to measure the associated coupled force and four isolated bio-potential preamplifiers (iwrox C-ISO-256) were put at the position of corresponding muscles as shown in Figure 3 to collect sEMG signals. As the two devices were separate, the NI USB-6210 device with 15,625 Hz per channel was used to make the gathered data simultaneous. During this experiment, the collected data-set contained the swing angle, velocity, acceleration data, the coupled force information in X-direction, and the real-time sEMG signals of the subject’s lower limb joints.

The experiment was also divided into two steps, the first for thigh swing data collection and the second for shank. Before each step, measuring the height from belt to joint position, that is, the length of the force arm as Figure 4, the joint toque can therefore be calculated as follows:

$${\overline{\mathbf{\tau }}}_{\mathrm{j}}={\overline{\mathbf{F}}}_{\mathrm{j}}{h}_{\mathrm{j}},$$

(11)

where ${\overline{\mathbf{\tau }}}_{\mathrm{j}}$ , ${\overline{\mathrm{F}}}_{\mathrm{j}}$ , and $h_{\mathrm{j}}$ is the calculated joint torque, measured force, and the length of the force arm in the experimental first step (j = t) and second step (j = s). As long as $h_{\mathrm{j}}$ has been determined, the subject performed swing motion using joints power at different initial angles by securing thigh of the human body in the first step and shank in the second to the experimental platform through the belt, as Figure 5(a) and (b). During the experiment, $h_{\mathrm{j}}$ remains constant, that is to say, the subject’s thigh and shank always keeps perpendicular to the ground separately in the first and second step. Notably, the movements are always within the threshold of this stiff platform. By comparing the predicted joint torque ${\hat{\mathbf{\tau }}}_{\mathrm{j}}$ and the joint torque ${\overline{\mathbf{\tau }}}_{\mathrm{j}}$ , the BPNN can therefore be evaluated.

Figure 4.

The measurements schematic of force arm length of: (a) thigh $h_{\mathrm{t}}$ and (b) shank $h_{\mathrm{s}}$ .

Figure 5.

The experimental movements with different initial angles of collecting (a) thigh and (b) shank swing data for network training at the stiff platform.

The trained neural network prediction

In order to validate the performance of the trained BPNN to predict joint torque from sEMG signals, a new data set should be utilized. Four Trigno Avanti sensors were used to collect sEMG signals and kinematic data simultaneously in Delsys system ⁴² while one put on the shank just for kinematic data collection with a frequency of 2000 Hz, which was shown in Figure 1(c). Freely swing motions were made as Figure 6 by the subject without exoskeleton throughout validation data collection. Information such as acceleration and rotation relayed by the sensors can be leveraged to discern movement activity time-synchronized with the sEMG signals, resulting in no time skew between data.

Figure 6.

The whole experiment movements in different time of collecting thigh and shank swing data simultaneously for network prediction at Delsys system.

Since the system can only collect velocity information, it needs to be differentiated to get the acceleration and integrated to acquire the angle. From the measured kinematic data $\overline{\mathit{\theta }}$ , the corresponding human joint torque during this swing experiment ${\overline{\mathbf{\tau }}}_{\mathrm{h}}$ can be computed using inverse dynamic analysis (IDA) according to equation (1) $({\overline{\mathbf{\tau }}}_{\mathrm{EXT}}=0)$ . After giving the feature of sEMG signals measured in this experiment into the trained BPNN, the predicted joint torque ${\hat{\mathbf{\tau }}}_{\mathrm{h}}$ will export. By comparing the two joint torque, the performance of the trained BPNN to predict joint torque from sEMG signals can be validate.

Data processing and analysis

Since the training subset data was captured at 15,625 Hz while the prediction subset data were sampled at 2000 Hz, the data were sub-sampled to their least common factor, 125 Hz, using a sliding window. ⁴⁰ There were four low-pass filters designed here to process measured data. The first is a finite impulse response filter with passband-frequency of 0.03 $\pi$  rad/s, stopband-frequency of 0.1 $\pi$  rad/s, passband-ripple of 0.5 dB, and stopband-attenuation of 65 dB to filter all data in Experiment (a). The second is a zero-lag fourth order butterworth filter to filter the force data in Experiment (b) at 10 Hz. The third is a butterworth filter with a passband-ripple of 1 and a stopband-attenuation of 30 to filter the sEMG signals in Experiment (a)/(b). The forth is a zero-lag fourth-order butterworth filter with a passband cutoff frequency of 6 Hz to filter the kinematic data in Experiment (c).

Three metrics were applied here to evaluate performance of the results, which contains the mean absolute error (MAE), the root-mean-square error (RMSE), and the Pearson correlation coefficient (PCC):

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{MAE}=\frac{1}{n}\sum _{\mathrm{s}=1}^{\mathrm{n}}|y_{\mathrm{s}}-{\hat{y}}_{\mathrm{s}}|,}\\ \multicolumn{1}{c}{\mathrm{RMSE}=\sqrt{\frac{1}{n}\sum _{\mathrm{s}=1}^{\mathrm{n}}{(y_{\mathrm{s}}-{\hat{y}}_{\mathrm{s}})}^2},}\\ \multicolumn{1}{c}{\mathrm{PCC}=\frac{\underset{(y_{\mathrm{s}},{\hat{y}}_{\mathrm{s}})}{\mathrm{cov}}}{\sigma y_{\mathrm{s}}\sigma {\hat{y}}_{\mathrm{s}}},}\end{array}$$

(12)

where $y_s$ and ${\hat{y}}_s$ is the expected and calculated values respectively, $n$ is the samples number, $\mathrm{cov}(y_{\mathrm{s}},{\hat{y}}_{\mathrm{s}})$ is the covariance matrix between the expected and calculated values. $\sigma y_{\mathrm{s}}$ and $\sigma {\hat{y}}_{\mathrm{s}}$ is the standard deviation of $y_{\mathrm{s}}$ and ${\hat{y}}_{\mathrm{s}}$ respectively. The smaller value of MAE and RMSE, as well as the larger value of PCC mean the better performance of the experimental results.

Results

Identification results

First, the desired excitation trajectory ${\mathit{\theta }}_{\mathrm{e}}$ was generated using the PSO method mentioned before. According to (6), the fitness could be calculated, which was shown in Figure 7(a). The value of fitness function is decreasing with increasing number of iterations. Taking the parameters computed at the point when the fitness function value is minimum, Figure 7(b) displays the trajectories of the knee and hip joints in a period.

Figure 7.

The result of PSO: (a) the value of fitness function (blue line) and the minimum value adopted (red rectangle) and (b) the trajectories of the knee and hip joints in $t\in [0,2\pi ]s$ period.

After importing the computed trajectories, the experiment can then be carried out, where parameters of unloaded exoskeleton robot ${\mathrm{X}}_{\mathrm{r}}$ and exo-human ${\mathrm{X}}_{\mathrm{hr}}$ were identified sequentially. Using ${\mathrm{X}}_{\mathrm{hr}}-{\mathrm{X}}_{\mathrm{r}}$ , the subject’s parameters ${\mathrm{X}}_{\mathrm{h}}$ can be computed, which are listed in Table 3. It can be seen that the identification results calculated is close to the average parameters of the human body. ²⁷

Table 3.

The experimental result and the statistically average data of identification parameters.

Parameters	Robot	Exo-human	Subject	Average
$X_1$	5.1946	9.0303	3.8357	2.0575
$X_2$	1.2842	3.7551	2.4709	0.6808
$X_3$	0.572	2.4405	1.8685	1.5761
$X_4$	2.4883	5.636	3.1477	3.8228
$K_{S1}$	3.7575	3.7575	0	0
$K_{m1}$	4.5604	4.5604	0	0
$K_{S2}$	4.1347	4.1347	0	0
$K_{m2}$	8.8975	8.8975	0	0

The actual joint torque ${\overline{\mathbf{\tau }}}_{r1}/{\overline{\mathbf{\tau }}}_{r2}$ and the computed joint torque ${\hat{\mathbf{\tau }}}_{r1}/{\hat{\mathbf{\tau }}}_{r2}$ mentioned before are compared in Figure 8. There is little difference in the trend between the computed and real torque levels. The absolute error between the computed and actual joint torques are also shown here in Figure 9. It can be seen that the error in human knee torque is the largest, which can up to about −40, however the others are within $\pm 20$ .

Figure 8.

The joint torque in the identification experiments: (a) the exoskeleton hip torque, (b) the exoskeleton knee torque, (c) the subject hip torque, and (d) the subject hip torque.

Figure 9.

The absolute joint torque error in the identification experiments of (a) exoskeleton robot and (b) subject.

Moreover, three metrics in equation (12) are listed in Table 4. The MAE and RMSE is low, with the largest RMSE value of 10.2 between human knee calculated and experimental joint torque. The PCC is high, with the lowest value of 0.84. Generally speaking, the identification performance of the exoskeleton robot is superior to the subject and the hip joint is better than the knee. Though the subjective body movements may produce influence, the accuracy and reliability of the identified parameters can still be validated from the above results. There exists a solid sign of identification combining the proximity of the result to the average parameters of human body and little error criteria between the actual measured joint torque to the computed value. Therefore, the identified parameters can be used to calculate the joint torques is assumed.

Table 4.

Three metrics in identification results.

Joint torques	MAE	RMSE	PCC
Exoskeleton hip joint	1.0757	1.5546	0.9842
Exoskeleton knee joint	2.6393	3.9527	0.8511
Human hip joint	4.8840	6.3289	0.9484
Human knee joint	6.9535	10.2334	0.8413

BPNN training results

According to equation (10), the hidden layer number is determined as 42 of the hip joint BPNN and 45 of the knee joint BPNN. Then the two constructed BPNN is trained until the training error cannot be reduced for six consecutive training sessions. When terminated, the actual iteration number, highest MSE, lowest MSE, gradient, and $\mu$ is 112, 3.02, 0.00259, 0.0835, $1\times {10}^{-5}$ for hip and 114, 4.52, 0.000625, 0.0122, $1\times {10}^{-5}$ for knee BPNN respectively.

Figure 10(a) and (b) illustrates 100 random sampling points of the output datasets and the corresponding prediction results of the network. It can be seen from Figure 10(c) and (d) that the error between the actual measured and the predicted joint torque is small. The error of hip BPNN is almost within $\pm 15$ and the error of knee BPNN is no more than $\pm 10$ . Besides, the visualization of regression capabilities of the trained BPNN are shown in Figure 11, the regression coefficient are both over 0.9. From the results all above, the two BPNN have been trained in a desirable performance.

Figure 10.

The joint torque measured actually and predicted by BPNN of (a) hip and (b) knee. The absolute error between the measured and predicted joint torque of (c) hip and (d) knee.

Figure 11.

The visualization of the regression capabilities of (a) hip and (b) knee BPNN.

BPNN prediction results

In the estimation trials, the inverse dynamics analysis (IDA) could be carried out through equation (4) using the testing set’s kinematic data $\overline{\mathit{\theta }}$ , which are shown in Figure 12(a) and (b) as well as the human body parameters ${\hat{\mathrm{X}}}_{\mathrm{h}}$ gleaned from the identification trials. According to a comparison between the joint torques computed by the IDA and predicted by the trained BPNN, the fitting impact of the trained model could be determined. The comparison results are illustrated in Figure 12(c) and (d). It is evident from the figure that the computed and estimated joint torques almost overlap, demonstrating the ability of neural networks with proper training to determine a subject’s intention for active motion in real time. To assess how well the estimated joint torque matched the theoretical real value, the absolute error between them are shown in Figure 12(e) and (f) and it can be seen that the error are almost within $\pm 15$ both for hip and knee joint.

Figure 12.

The measured joint angular velocity of (a) hip and (b) knee. The joint torques computed by IDA and estimated by BPNN of (c) hip and (d) knee. The absolute error between the computed and estimated joint torque of (e) hip and (f) knee.

In addition, three calculated metrics through (12) are listed in Table 5. The MAE and RMSE between IDA computed and BPNN predicted joint torque are all within 5, where the largest is the RMSE in hip joint estimation with the value of 4.93. The small value of MAE and RMSE indicate the trained model have a desirable performance. However, the PCC are not so high because the static joint torque measured in Experiment (b) may not so relevant with the dynamic joint torque calculated using data from Experiment (c). From the validation results above, the effectiveness of trained BPNN can be proved, the correctness of identification parameters can also be counter-validated.

Table 5.

Three metrics in estimation results.

Joint torques	MAE	RMSE	PCC
Hip joint	2.8593	4.9307	0.7908
Knee joint	1.4236	3.9327	0.6510

Discussion

To ensure effective and safe human-exoskeleton interaction – a critical component of intelligent rehabilitation – as well as to enable robots to assist subjects in human-exoskeleton collaboration, accurate and robust joint torque estimation during human movement is necessary, which can also be utilized to control rehabilitation exoskeleton robots during complex movements appropriately. Moreover, the output estimation method can expand into fields of humanoid robots and even any continuous-time systems while not limited to exoskeleton robots, as determining the control architecture for them all relies on the output-feedback.^43,44 Surface electromyography (sEMG) signals have been shown to be a reliable method of assessing human intended motions in coupled human-exoskeleton systems. However, there exist two problems in using sEMG signals to estimate human joint torque. First, the subject plays a critical role in human-exoskeleton coupling exercises, which means the average human body parameters used in human-exo coupling dynamics modeling are not exact. Second, the conventional techniques including creating biomechanical or regression models that estimate continuous joint motion variables are too sensitive to changes in the input signal. That is to say, the instability property of sEMG signals makes it difficult to reliably predict continuous joint motion, significantly restricting the method’s generalizability. This study proposes a complete process contains three main parts as Figure 1 for human dynamics analysis based on solving the above two problems.

In this paper, we performed human parameter identification experiments prior to human joint torque estimation for the first time. Following the identification of the human body parameters, we used the human lower limb joint torque dataset and the sEMG recorded during the swing experiment to train the BP neural network we proposed. In the end, the human lower limb joint torque estimation was performed using the trained BP neural network model. The outcomes demonstrate that this procedure is able to precisely identify human parameters and provides an acceptable approximation of the joint torques involved in the lower limb swing of an individual.

In preparation for the identification trials, we first created the human-exoskeleton coupled dynamics model, as seen in Figure 2. The most suitable ideal trajectory for swing in the hip and knee joints was determined by applying the PSO algorithm satisfying (6), and the results are displayed in Figure 7(b). Then the unloaded exoskeleton robot platform used this ideal trajectory as input for the swinging experiment. The exoskeleton robot’s parameters could be identified using both the actual joint trajectories and the driving torques that were gathered in the swing experiment. The human body must be propelled by the exoskeleton robot during the swinging experiment in order to evaluate the parameters of the human body. We utilized the final identified data, which are displayed in Table 3 to compute joint torque in equation (7), and the comparison of it and the measured joint torque is seen in Figure 8. The validity of our identification studies may be demonstrated by the similarity of the calculated and real torques as well as the identified data and the average human body data.

For the BPNN training process, the data-set is made up of the synchronized sEMG produced after fixating the human lower limb with the human joint torque data computed by equation (9). Using the Neural Network Toolbox of MATLAB, the results are shown in Figure 10. For the BPNN validation process, the data-set contains the sEMG signals, which were fed into the trained BP neural network and kinematic data of the human lower limb swinging freely, which were used with the identified parameters to calculated the corresponding joint torques. The comparison between the joint torques from BPNN estimated and those derived from IDA calculated, as revealed by Figure 12, was evident. The accuracy the neural network model’s prediction was could be confirmed by the small mean absolute error (MAE) and root-mean-square error (RMSE) in Table 5, which also show reversely the correctness of the identification parameters.

A limitation of this study is that we only depict the experimental process, another influencing factors such as the control algorithm of exoskeleton robot is not considered. The impedance control method may become our future research directions, as it is superior in the robotic systems with contact force.^45,46 Another limitations are about the experiments, for the parameter identification process (Experiment (a)), the human body should be fully driven by the driving moments generated by the exoskeleton robot to swing, but the subject in this experiment only subjectively did not exert any force, we nevertheless did not use electromyography-types of equipment for confirmation. A further limitation of this study is the placement of equipment for measuring sEMG signals, which we were unable to analyze anatomically and had to rely on visual observation of muscle orientation to determine. Furthermore, the limits of the experimental apparatus precluded us from validating the BP neural network’s dynamic torque estimate accuracy. Even though we have enough experimental data for the current investigation, it is limited to one participant and hasn’t been developed into a data-set group with individual variability. Besides, due to the limitations of our experimental equipment, the prediction quality is not so satisfying, which, however, we are struggling to improve.

Conclusion

In this paper, a complete process for analyzing dynamics of human motion using movements data and surface EMG sigals is presented. Human parameters for subject can be identified from human-exoskeleton coupling experiments, and human joint torques can be estimated by a trained BP neural network. The process is comprehensive and the results can be cross-checked at each step. The proposed identification method can be verified from little errors between the identified and averaged human body parameters, and the measured and computed joint torques. Moreover, the accuracy of the prediction of BPNN can not only demonstrate the successful of network training but also the correctness of identification results. The above-described analysis procedure can be extended to various types of human motion movement.

In further work, we will examine the application of the proposed protocol in various subjects and assess the suggested approach in more practical contexts, such as walking on flat ground and climbing stairs. We will also do our best to improve our experimental equipment and the variegation of control methods.