Gaussian Process Autoregression for Joint Angle Prediction Based on sEMG Signals

There is uncertainty in the neuromusculoskeletal system, and deterministic models cannot describe this significant presence of uncertainty, affecting the accuracy of model predictions. In this paper, a knee joint angle prediction model based on surface electromyography (sEMG) signals is proposed. To address the instability of EMG signals and the uncertainty of the neuromusculoskeletal system, a non-parametric probabilistic model is developed using a Gaussian process model combined with the physiological properties of muscle activation. Since the neuromusculoskeletal system is a dynamic system, the Gaussian process model is further combined with a non-linear autoregressive with eXogenous inputs (NARX) model to create a Gaussian process autoregression model. In this paper, the normalized root mean square error (NRMSE) and the correlation coefficient (CC) are compared between the joint angle prediction results of the Gaussian process autoregressive model prediction and the actual joint angle under three test scenarios: speed-dependent, multi-speed and speed-independent. The mean of NRMSE and the mean of CC for all test scenarios in the healthy subjects dataset and the hemiplegic patients dataset outperform the results of the Gaussian process model, with significant differences (p < 0.05 and p < 0.05, p < 0.05 and p < 0.05). From the perspective of uncertainty, a non-parametric probabilistic model for joint angle prediction is established by using Gaussian process autoregressive model to achieve accurate prediction of human movement.

1. Introduction

Rehabilitation robots and other rehabilitation equipment have developed rapidly and are widely used for therapeutic training of patients suffering from neurological disorders including stroke, cerebral palsy and spinal cord injury. Patients with impaired neurological function are able to use rehabilitation equipment for a variety of exercises to restore strength and flexibility in their extremities (1). Neurorehabilitation techniques can help both passive and active training of patients with neurological injuries. Compared to passive training through rehabilitation equipment, training that involves the patient's own will can improve the effectiveness of treatment and restore motor function through active movement (2). Electromyography (EMG) represents the sum of subcutaneous motor action potentials generated through muscle contraction (3), which can represent neuromuscular activity and is a way to reflect the patient's voluntary effort. Surface electromyography (sEMG) is a non-invasive EMG signal that is mainly acquired through non-invasive electrodes, which has the advantage of being more accessible and less likely to impede the normal activity of the user and will override the actual joint motion (4), and is often used as a control command to realize the human-robot interface (HRI) of portable or wearable assisted rehabilitation equipment. The surface electromyography (sEMG) is a non-invasive EMG signal collected by surface electrodes, which has the advantage of being more accessible than electroencephalography (EEG) and less likely to impede the user's normal activities. The forces generated by muscles in response to neural control signals depend on a large number of variables distributed over many spatiotemporal scales (5), which makes it difficult to predict the muscle force. While the nervous system clearly has knowledge of some of these variables, such as muscle length and velocity, other aspects of muscle dynamics (e.g., the detailed dynamics of molecules at the level of individual myofibrils and sarcomeres) are much more difficult to measure and estimate. With knowledge of the nervous system, it is possible to model the relationship between neural control signals and muscle force and use it to predict or simulate muscle force production. However, this relationship itself cannot be accurately estimated, and the remaining unaccounted aspects of muscle dynamics will result in seemingly random fluctuations in force, also known as motor noise. Therefore, there are two main sources of uncertainty in the neuromusculoskeletal system: irreducible noise during the motor system and variability in the relationship between neural control signals and muscle outputs (6). The impact of uncertainty in neuromusculoskeletal models on joint motion prediction can be mitigated, but not completely eliminated, by different modeling approaches. The human neuromuscular system is a highly non-linear and time-varying system with uncertainty (7). In order to allow models to handle the dynamic high-dimensional nature of the neuromuscular system, it is not enough to rely on traditional model structures and faster computational processing. Therefore, a central issue for further research on neuromusculoskeletal systems, or any artificial controller, is how to command muscles effectively in the presence of uncertainty.

Non-linear system modeling and identification can be divided into parametric and non-parametric models from the perspective of model structure (8–12). From the perspective of the Bayesian statistical framework, probability distributions over the functional space can be considered and modeled by optimizing these distributions to characterize uncertainty. This type of model has no explicit modeling mechanism or constraints and is referred as non-parametric modeling (13). The non-parametric model does not depend on a specified set of parameters or a fixed model structure, and is an estimation method based on statistical principles, which can generate functions to fit the data without the constraints of the model mechanism. The number and nature of the parameters of the non-parametric model are flexible and variable, which can change accordingly with the change of the data set. It can be used for modeling and analyzing high-dimensional time-varying systems, taking the uncertainty into account during the modeling process, and characterizing the uncertainty. The non-parametric methods mainly include spectral estimation, spectral analysis, correlation analysis and kernel-based analysis methods. Gaussian process (GP) model is a non-parametric kernel method in the framework of Bayesian model, which is simple to implement, computationally efficient, and most importantly, GP model can describe the posterior distribution of the model function, which in turn can be used to describe the uncertainty of the model, and is a recent research hotspot for non-parametric methods (14, 15). Kang et al. proposed an effective method for generating suboptimal motion of a multi-body system using a GP dynamics model to achieve dimensionality reduction of the system and deal with motion optimization problems (16). Schearer and Ullauri achieved the estimation of joint moments by building semi-parametric and non-parametric models through GP model, respectively (17, 18). Xiloyannis et al. used a Gaussian autoregressive model for decoding neural information to achieve multidimensional decoding control of 11 joint movements (19). Yang et al. proposed a proportional pattern recognition control of arm muscles using a wearable ultrasound sensor to achieve both gesture recognition and muscle contraction force estimation based on statistical features and Gaussian process regression models (20).

It has been shown that joint motion can reflect the intrinsic dynamics of human movement (21), so motion signals can also be used in the modeling of the neuromusculoskeletal system for building autoregressive models for prediction. The non-linear autoregressive with eXogenous inputs (NARX) model is an effective method for solving non-linear sequential problems, and modeling in conjunction with the NARX model can better incorporate the non-linear spatiotemporal correlation structure of muscle-driven control signals and natural human motion (22). Dynamic recurrent neural networks based on the NARX model are widely used for joint angle estimation, decoding shoulder, elbow and wrist motions and prosthesis model control (23–25). Gupta proposed an ankle joint angle estimation model based on the NARX model using sEMG signals and knee joint angle signals, and the performance of the model proved its applicability to ankle joint angle estimation for active prosthesis, orthosis and lower limb rehabilitation (26). Liu et al. used the NARX model to train and identify the EMG signals motion mapping relationship between a rehabilitation training bed and sEMG signals based motion prediction to achieve the identification of upper body tilt in different directions (27). Raj et al. proposed a multilayer perceptron neural network model based on the NARX model for estimating elbow joint angle and elbow joint angular velocity, and the proposed model estimated elbow joint angle and elbow joint angular velocity with high accuracy (28).

Since the neuromusculoskeletal system is a dynamic time-varying system, the GP model is only a mapping of the input to the output distribution, which is a basic static system. By combining the NARX model into the modeling, the resulting model can not only adapt to non-linear discrete-time processes, but also to different noise models, and the resulting dynamic model, which will better fit the physiological properties of the neuromusculoskeletal system. Thus, in this paper, a knee joint angle prediction model based on sEMG signals is proposed by considering the physiological properties of microscopic muscle activation and combining NARX with GP model. First, considering the physiological properties of sEMG, the muscle activation kinetic model is used to extract features from the sEMG signals of a pair of antagonist muscles controlling knee joint motion, and the muscle activation intensity of this pair is obtained. Then a joint angle prediction model based on GP model is proposed for the uncertainty of the neuromusculoskeletal system with the muscle activation intensity as the input signal. Since the neuromusculoskeletal system is a time-varying non-linear system, the one-step ahead (OSA) prediction of the NARX model is introduced to construct a Gaussian process autoregressive model, which uses the confidence interval of the prediction to describe the uncertainty, reduces the influence of model uncertainty on the prediction results, and improves the prediction rationality, accuracy, and efficiency of the joint angle prediction model.

2. Materials and Methods

2.1. Muscle Activation Dynamics

The sEMG signal is the sum of action potentials recruited to the muscle by surface electrodes and is used to reflect the activation level of the muscle. sEMG signals can be considered as a form of characterization of neuromotor control commands and are widely used to analyze musculoskeletal models. The feature extraction of the EMG signals by the existing studied models only considers the macroscopic characteristics of the EMG signals, without considering the microscopic physiological characteristics of muscle activation. In order to characterize the time-varying features of the sEMG signals, respond to micro-physiological properties, and reflect the relationship between EMG signals, neural activation and muscle activation, a muscle activation kinetic model was established to achieve feature extraction of the sEMG signals (29–31). The muscle activation kinetics is mainly expressed as the transformation process between EMG signals and muscle activation, as shown in Figure 1, where a(t) is the muscle activation, e(t) is the processed sEMG signal, q(t) the neural activation, and detailed in Li et al. (32).

FIGURE 1

Figure 1. Muscle activation dynamics.

2.2. Gaussian Process

Gaussian process (GP) is defined as a random process consisting of infinite high-dimensional random variables in a high-dimensional space, in which the joint distribution among any finite number of random variables is a Gaussian distribution. GP model can be derived from the weight-space view or the function-space view. Since each set of weights implies a specific function and the distribution of the weights implies the distribution of the function, the distribution of the GP can be obtained from the function-space view to obtain the equivalent of the weight-space view (33), which is the more commonly used derivation method for Gaussian process models.

Suppose the sample set ${D}$ has N samples:

$\begin{array}{ll}{D}=(\text{X},\text{Y})=\left\{({\text{x}}_i,y_i)|{\text{x}}_i\in {ℝ}^d,y_i\in ℝ,i=1,\dots ,N\right\}& (1)\end{array}$

where x_i denotes input vector,y_i denotes output vector.

A Gaussian process is completely specified by its mean function m(x) and covariance function k(x, x′) (34):

$\begin{array}{ll}f(\text{x})~{G}{P}(m(\text{x}),k(\text{x},{\text{x}}^{\prime }))& (2)\end{array}$

$\begin{array}{ll}m(\text{x})=𝔼\left[f(\text{x})\right]& (3)\end{array}$

$\begin{array}{ll}k(\text{x},{\text{x}}^{\prime })=𝔼\left[(f(\text{x})-m(\text{x}))(f({\text{x}}^{\prime })-m(\text{x}))\right]& (4)\end{array}$

where the random variable function f(x) represents the distribution of y_i at x_i, the mean function m(x) reflects the expected function value at input x. The covariance function models the dependence between the function values at different input points x and x′, which is often referred as the kernel function of a GP model.

In the Gaussian process regression, considering the following model:

$\begin{array}{ll}\text{Y}=f(\text{X})+\epsilon & (5)\end{array}$

where X denotes the input vector,Y denotes the observed vector with noise, the noise follows a Gaussian distribution $\epsilon ~{N}(0,{\sigma }^2)$, the random variable function f(X) follows a Gaussian distribution:

$\begin{array}{ll}f(\text{X})~{N}(\mu (\text{X}),\text{K}(\text{X},\text{X}))& (6)\end{array}$

Thus:

$\begin{array}{ll}\text{Y}~{N}(\mu (\text{x}),\text{K}(\text{X},\text{X})+{\sigma }^{2}I)& (7)\end{array}$

For the prediction input $X^{\text{*}} \text{=} {(x_1^{*},\dots ,x_N^{*})}^T$, the joint distribution of the predicted values f(X^*) and the training data output is:

$\begin{array}{ll}\left(\begin{array}{c}\text{Y}\\ f({\text{X}}^{*})\end{array}\right)~{N}\left(\left[\begin{array}{c}\mu (\text{X})\\ \mu ({\text{X}}^{*})\end{array}\right],\left[\begin{array}{cc}\text{K}(\text{X},\text{X})+{\sigma }^2\text{I}& \text{K}(\text{X},{\text{X}}^{*})\\ \text{K}({\text{X}}^{*},\text{X})& \text{K}({\text{X}}^{*},{\text{X}}^{*})\end{array}\right]\right)& (8)\end{array}$

where K(X, X) denotes the covariance matrix of the input signal, which is a symmetric semi-positive definite matrix of N × N order:

$\begin{array}{ll}\text{K}(\text{X},\text{X})=\left(\begin{array}{ccc}k({\text{x}}_1,{\text{x}}_1)& \cdots & k({\text{x}}_1,{\text{x}}_N)\\ ⋮& \ddots & ⋮\\ k({\text{x}}_N,{\text{x}}_1)& \cdots & k({\text{x}}_N,{\text{x}}_N)\end{array}\right)& (9)\end{array}$

Knowing the joint high-dimensional distribution, the posterior distribution is obtained by finding the conditional probability p(f(X^*)|Y, X, X^*) according to Bayes' theorem:

$\begin{array}{ll}f({\text{X}}^{*})|{\text{X}}^{*}~{N}({\mu }^{*},{\Sigma }^{*})& (10)\end{array}$

$\begin{array}{ll}{\mu }^{*} \text{=} K(X^{*},X)[k{(X,X+{\sigma }^{2}I]}^{-1}(Y-\mu (X))+\mu (X^{*})& (11)\end{array}$

$\begin{array}{ll}{\Sigma }^{*}=\text{K}({\text{X}}^{*},{\text{X}}^{*})-\text{K}({\text{X}}^{*},\text{X}){[\text{K}(\text{X},\text{X})+{\sigma }^2\text{I}]}^{-1}\text{K}(\text{X},{\text{X}}^{*})& (12)\end{array}$

The covariance matrix accounts for the major part of the posterior distribution of the Gaussian process, and the covariance matrix is a key component of the Gaussian process prediction. The kernel function is the main structure of the covariance matrix, so it is also a central part of the Gaussian process model. In the modeling and identification of dynamic systems, the dimensions of the inputs are relatively high, which makes the description of the mapping function complex. Considering the smoothness and continuity of dynamic systems, the squared exponential (SE) kernel function is often used in the modeling process (35), which is defined as:

$\begin{array}{ll}k(\text{x},{\text{x}}^{\prime })={\sigma }_f^{2}exp\left[-\frac{{(\text{x}-{\text{x}}^{\prime })}^T(\text{x}-{\text{x}}^{\prime })}{2{\sigma }_l^2}\right]& (13)\end{array}$

where hyperparameter σ_l is the characteristic length scale, which determines the relative weights of the distances of the input variables. σ_f is the signal standard deviation, which reflects the magnitude of the function change.

The Gaussian process model is mainly determined by the kernel function and its hyperparameters, and its learning process is a process of training through the data to obtain the posterior probability distribution, which mainly includes the selection (or design) of the kernel function and the determination of the hyperparameters.

The kernel functions of Gaussian processes often contain unknown and indefinite hyperparameters, such as length scales, signal and noise variances, etc. These need to be inferred from the data, resulting in posterior distributions of the hyperparameters that are not easily obtained. Therefore, the full Bayesian derivation of hyperparameters is not commonly used in practical applications. The usual practice is to obtain point estimates of the hyperparameters by maximizing the log marginal likelihood.

Given a sample set ${D}$ and the hyperparameter of the Gaussian process is θ, the marginal likelihood is as shown in Equation (14):

$\begin{array}{ll}p(\text{Y}|\text{X},\theta )={\displaystyle \int }p(\text{Y}|\text{X},f,\theta )p(f|\text{X},\theta )df& (14)\end{array}$

The marginal likelihood is mainly a marginalization of the function. In the Gaussian process model, the prior f|X, θ of the model is a Gaussian distribution, i.e., $p(f|\text{X},\theta )={N}(0,{\text{K}}_{\theta })$. When the observed likelihood function p(Y|X, f, θ) of the sample set is also Gaussian distributed, i.e., $p(\text{Y}|\text{X},f,\theta )={N}(f,{\sigma }^2\text{I})$, then p(Y|X, θ) is also Gaussian distributed:

$\begin{array}{ll}p(\text{Y}|\text{X},\theta )={\displaystyle \int }{N}(0,{\text{K}}_{\text{Y}}){N}(f,{\sigma }^2\text{I})df={N}(0,{\text{K}}_{\text{Y}}+{\sigma }^2\text{I})& (15)\end{array}$

According to Equation (15), the log marginal likelihood is obtained as:

$\begin{array}{ll}logp(\text{Y}|\text{X},\theta )=-\frac{1}{2}{\text{Y}}^T{{\text{K}}_{\text{Y}}}^{-1}\text{Y}-\frac{1}{2}log|{\text{K}}_{\text{Y}}|-\frac{N}{2}log2\pi & (16)\end{array}$

where ${\text{K}}_{\text{Y}}=\text{K}(\text{X},\text{X})+{\sigma }^2\text{I}$ is the output covariance matrix.

The maximum likelihood estimation combined with the conjugate gradient method is commonly used for the Gaussian process model to achieve the estimation of the hyperparameters of the model, and the computational complexity of this method is ${O}(N^2)$ for each hyperparameter, and the computational complexity is small. The hyperparameter estimates of the Gaussian process model are obtained by maximizing the log marginal likelihood function through a gradient ascent based optimization tool:

$\begin{array}{ll}\begin{array}{c}\frac{\partial }{\partial {\theta }_i}logp(\text{Y}|\text{X},\theta )=-\frac{1}{2}{\text{Y}}^T{{\text{K}}_{\text{Y}}}^{-1}\frac{\partial {\text{K}}_{\text{Y}}}{\partial {\theta }_i}{{\text{K}}_{\text{Y}}}^{-1}\text{Y}-\frac{1}{2}tr({{\text{K}}_{\text{Y}}}^{-1}\frac{\partial {\text{K}}_{\text{Y}}}{\partial {\theta }_i})\\ =\frac{1}{2}tr\left(({\alpha \alpha }^T-{{\text{K}}_{\text{Y}}}^{-1})\frac{\partial {\text{K}}_{\text{Y}}}{\partial {\theta }_i}\right)\end{array}& (17)\end{array}$

where $\alpha ={\text{K}}_{\text{Y}}^{-1}\text{Y}$.

2.3. Non-parametric Model for Joint Angle Prediction Based on sEMG Signals

2.3.1. Joint Angle Prediction Based on Gaussian Process Model

The hamstrings and quadriceps are antagonistic muscles that together control the flexion and extension of the knee joint. The hamstrings are the muscles of the posterior thigh and consist mainly of the semitendinosus, semimembranosus and biceps femoris, while the quadriceps are the muscles of the anterior thigh and consist mainly of the vastus lateralis, vastus medialis, vastus intermedius and rectus femoris. In this paper, the sEMG signals of a pair of muscles in this antagonistic muscle group, the semimembranosus and the lateral femoris, were selected for the development of a non-parametric model for joint angle prediction. The physiological properties of muscle activation are combined with the GP model, and the squared exponential kernel is selected for subsequent modeling and analysis to establish a non-parametric model for joint angle prediction based on the GP model, as shown in Figure 2, where k denotes kth time step, e_{1, k} and e_{2, k} are the preprocessed sEMG signals of semimembranosus and lateral femoris muscles, respectively, which are then subjected to muscle activation dynamics to calculate the muscle activation a_{1, k} and a_{2, k}. u_k denotes the input of the GP model, $u_k={[a_{1,k},a_{2,k}]}^T$. The output of the Gaussian process model ŷ_k is the predicted value of the joint angle.

FIGURE 2

Figure 2. Joint angle prediction based on Gaussian process model.

2.3.2. Joint Angle Prediction Based on Gaussian Process Autoregressive Model

Since the neuromusculoskeletal system is a dynamic time-varying system, the joint angle prediction based on GP model only maps the input to the output distribution, which is a basic static system. Therefore, this paper further considers the time series of input and output signals to establish dynamic Gaussian process regression.

Non-linear autoregressive with eXogenous inputs (NARX) model is an effective method for solving non-linear sequence problems that accommodates non-linear discrete time processes and noisy models (36, 37). NARX is a dynamic recurrent network that predicts the current value of the system output by using a non-linear function f with previous inputs and outputs. The NARX model based on the actual measured values of the output is called One-step ahead (OSA) prediction:

$\begin{array}{ll}{y}_k^{*}=f(u_k,u_{k-1}\dots ,u_{k-n_u},y_{k-1},y_{k-2},\dots ,y_{k-n_y})& (18)\end{array}$

where f(·) is the non-linear function between the input u_k and the estimated value $y_k^{*}$ and y_k denotes measured value of model output, k represents kth time step, n_u and n_y are the maximum lags for model input and output, respectively.

The joint angle signal can be easily collected by inertial measurement unit (IMU), etc., and the joint angle prediction system can be established by EMG signals to achieve further advance prediction of joint angle. Since the high accuracy of joint angle prediction is required in practical applications, this paper improves the joint angle prediction method based on GP model by using the NARX model and muscle activation dynamics, which establishes a Gaussian dynamic model with NARX structure, i.e., Gaussian process autoregressive model, for joint angle OSA prediction, as shown in Figure 3, where y_k denotes measured value of joint angle, ŷ_k denotes the joint angle prediction of Gaussian process autoregressive model, N_y and N_y are the maximum lags for model input and output, respectively.

FIGURE 3

Figure 3. Joint angle prediction based on Gaussian process autoregressive model.

3. Experiments and Results

3.1. Datasets

In this paper, publicly available datasets (dataset 1 and dataset 2) were cited to validate and analyze the proposed knee joint angle prediction model.

3.1.1. Healthy Subjects Dataset

Dataset 1 contains gait data from 10 healthy subjects in running condition, and Dataset 1 is published in https://simtk.org/projects/nmbl_running. The dataset measured knee moment signals, EMG signals and motion data of 10 healthy subjects running on a treadmill at four speeds (2.0, 3.0, 4.0, and 5.0 m/s). Subjects were all male (age: 29 ± 5 years; height: 1.77 ± 0.04 m; weight: 70.9 ± 7.0 kg), all provided informed consent, and each subject was experienced in long-distance running, at least 50 km per week. Fifty-four reflective markers were placed on each subject, and the trajectory of the markers was recorded using eight Vicon MX40+ cameras with a data acquisition frequency of 100 Hz. Ground reaction forces and moments were acquired using a Bertec Corporation treadmill with a sampling frequency of 1,000 Hz. Motion data and ground reaction forces were preprocessed with 4th order zero-phase hysteresis Butterworth low-pass filtering (cutoff frequency 15 Hz) and critical damping low-pass filtering (cutoff frequency 15 Hz), respectively. The Delsys Bangoli system was used to collect EMG signals. A total of 11 muscles including the gluteus maximus, biceps femoris long head, medial femoris, lateral femoris, semimembranosus, and tibialis anterior muscles were collected. Hamner et al. gave a complete description of the dataset (38).

3.1.2. Hemiparetic Subject Dataset

Joint angle prediction is mainly applied to the development of rehabilitation equipment and rehabilitation training, so a gait dataset containing a male patient with high-functioning right hemiparesis, which is publicly available at https://simtk.org/projects/emgdrivenmodel, was also selected for further testing and analysis of the knee angle prediction model proposed in this paper. The subject was 79 years old, height 1.7 m, mass 80.5 kg, with a LE Fugl-Meyer motor assessment score of 32/34 and right-sided hemiparesis. All experimental procedures were approved by the University of Florida Health Sciences Center Institutional Review Board (IRB-01), and the subject signed a written informed consent prior to participation in the experiment. The dataset collected gait data from subjects walking on a split-belt instrumented treadmill (Bertec Corp., Columbus, OH) at five different speeds (0.4, 0.5, 0.6, 0.7, and 0.8 m/s), with over 50 gait cycles collected for each speed. Motion capture in the experiments was mainly performed by an optical motion capture system (Vicon Corp., Oxford, UK) and ground reaction force detection was measured using the treadmill with sampling frequency of 100 and 1,000 Hz, respectively. The EMG signal acquisition was performed using Motion Lab Systems with a sampling frequency of 1,000 Hz. Ground response and marker motion data was filtered using a fourth-order zero-phase lag Butterworth filter with a cutoff frequency of 7 divided by the gait period. EMG signal data were collected for 16 muscle groups of the lower extremity, including the anterior tibialis, semimembranosus, long head of the biceps femoris, medial femur, and lateral femur. A full description of this dataset is provided by Meyer et al. (39).

3.1.3. Pre-processing and Feature Extraction Results

The sEMG signals of dataset 1 and dataset 2 are preprocessed. Raw sEMG signals was firstly filtered using Butterworth zero phase shift bandpass filter (4th order, cutoff frequency 40 Hz) to eliminate low-frequency noise, then full-wave rectified and low-pass filtered (4th order, cutoff frequency = 3.5/step period), and finally sEMG data from each muscle were normalized to the maximum value over all trials to obtain the processed sEMG signals e(t).

Feature extraction was then further performed using the muscle activation kinetic model in section 2.1. The results of pre-processing and feature extraction for dataset 1 (subject 1, 2 m/s) for the lateral femoral and semimembranosus muscles are shown in Figures 4, 5.

FIGURE 4

Figure 4. Feature extraction of the lateral femoral.

FIGURE 5

Figure 5. Feature extraction of semimembranosus.

3.2. Data Allocation Strategy

To study the effect of different speeds on joint angle prediction, joint angle prediction results were analyzed for two datasets of healthy subjects and hemiplegic subject in three conditions: speed-dependent, multi-speed and speed-independent. The normalized root mean square error (NRMSE) and correlation coefficient (CC) were used as evaluation indicators of the prediction performance of the joint angle, as shown in Equations (19) and (20).

$\begin{array}{ll}NRMSE=\frac{1}{y_{real}^{max}}\sqrt{\frac{1}{num}{\sum }_{i=1}^{num}{(y_{est}-y_{real})}^2}& (19)\end{array}$

$\begin{array}{ll}CC=\frac{cov(y_{est},y_{real})}{{\sigma }_{y_{est}}·{\sigma }_{y_{real}}}& (20)\end{array}$

where num denotes the number of samples tested, y_est denotes the predicted value of joint angle, y_real denotes the actual value of joint angle, and $y_{real}^{max}$ denotes the maximum magnitude of the actual joint angle value.

One-way analysis of variance (ANOVA) was conducted to assess the statistical difference of estimation errors obtained by different models (40). The level of statistical significance was set to p < 0.05.

Dataset 1 contained 10 subjects, each speed containing five gait cycles, and three conditions were analyzed for each subject. speed-dependent took the first three cycles of each speed as a training set, and the last two cycles at the corresponding speed as test set; multi-speed took the first three cycles of each speed together as a training set, and the last two cycles of each speed separately as a test set; the speed-independent took the last two cycles of one speed in turn as test set, and the first three cycles of each unselected speed together as training set.

Dataset 2 contained the left and right leg gait data of a subject with right-sided hemiplegia, and the left and right leg gait data were analyzed for three conditions. speed-dependent tested the angle prediction results of each speed, using 10 gait cycles of a single speed as the training set and another 10 cycles of the same speed as the test set; multi-speed used data from 5 m/s, each with 10 gait cycles, for a total of 50 gait cycles as the training set, and 10 gait cycles for each speed as the test set; speed-independent test took 10 gait cycles of one speed in turn as the test set, and 10 gait cycles of each of the other unselected speeds together as the training set.

Before performing the joint angle prediction based on Gaussian process model and Gaussian process autoregressive model, the Gaussian process model needs to be trained offline for different datasets according to the data allocation strategy. In this paper, the dataset was trained and tested on the MATLAB platform, and the Gaussian process regression model was trained offline using the “fitrgp” function. The method of estimating the model parameters was set to “exact,” and the point estimates of the hyperparameters were obtained by maximizing the log marginal likelihood, and the kernel function was set to the squared exponential kernel function. The input signal for offline training of joint angle prediction based on Gaussian process model was the muscle activation of lateral femoral and semimembranosus muscles, and the output was the normalized joint angle signal. Given that muscle dynamics is a second-order model, the joint angle prediction based on the Gaussian process autoregressive model is set to second order, which lead to the number of the maximum lags for model input and output be 2, i.e., n_u = n_y = 2. Therefore, the output signal of offline training was the joint angle signal, and the input was the muscle activation. The model was trained offline and joint angle prediction was performed according to the data allocation strategy for different datasets under different data allocation strategies, respectively.

3.3. Estimation Results

3.3.1. Joint Angle Prediction Results for Dataset 1

The proposed method was tested using the data in dataset 1. Subjects were tested for knee angle prediction in three cases, speed-dependent, multi-speed, and speed-independent, according to the data allocation strategy, and the test results for subject 9 knee angle prediction were shown in Figures 6–8, where “NARX-GP” was the joint angle prediction based on the Gaussian process autoregressive model, “GP” was the joint angle prediction based on the Gaussian process model, and “measurement” was the actual measurement of joint angle. The gray shading was the 95% confidence interval (μ±2σ) for the prediction of the joint angle based on the Gaussian process autoregressive model to describe the uncertainty. From Figures 6–8, it can be concluded that the direct joint angle prediction by Gaussian process model cannot describe the relationship between sEMG signal and joint angle well, and the knee joint angle prediction results have a large error. Establishing a Gaussian process autoregressive model for OSA prediction of joint angle can significantly improve the prediction accuracy and can approximate the actual joint angle signal. OSA prediction incorporates the actual values of the previous moments of output into the model structure with high prediction accuracy, and is suitable for scenarios where the actual measurements of the output are easy to collect and where high prediction accuracy is required. The joint angle signal can be easily collected by inertial measurement unit (IMU), etc., and the joint angle prediction system can be established by EMG signals to achieve further advance prediction of joint angle.

FIGURE 6

Figure 6. Joint angle prediction under speed-dependent of subject 9.

FIGURE 7

Figure 7. Joint angle prediction under multi-speed of subject 9.

FIGURE 8

Figure 8. Joint angle prediction under speed-independent of subject 9.

Further error assessment and statistical analysis of the prediction results were performed, and the mean NRMSE and CC between the predicted and actual measurements for different velocity joint angles of the subjects are shown in Figures 9, 10. It can be seen from the figures that the prediction results of the NARX-GP model were significantly better than those of the GP model, the NRMSE between the prediction results of the NARX-GP model for knee joint angle and the actual knee joint angle was small and significantly smaller than that of the GP model, and the strong correlation between the prediction results of the NARX-GP model and the actual values of the joint angle with a higher correlation coefficient than that of the GP model.

FIGURE 9

Figure 9. Average NRMSE for prediction of joint angle under different speed.

FIGURE 10

Figure 10. Average CC for prediction of joint angle under different speed.

The means and standard deviations of NRMSE and CC between predicted and actual values of joint angles for all subjects in speed-dependent, multi-speed and speed-independent conditions are shown in Tables 1, 2. From Tables 1, 2, it can be seen that the predictions of the NARX-GP model were highly correlated, and the NRMSE of the predictions was significantly lower than that of the GP model. The mean NRMSE was further calculated as 0.0039 ± 0.019, 0.0038 ± 0.0019, and 0.0059 ± 0.0061 for all subjects in the NARX-GP model in the load-dependent, multi-load and load-independent conditions, 0.1728 ± 0.0543, 0.1795 ± 0.0563, and 0.2003 ± 0.0659 for the GP model. ANOVA was used to evaluate NRMSE of NARX-GP and GP model predictions, showing significant differences (p < 0.05, p < 0.05, and p < 0.05). The NARX-GP model prediction errors were smallest for all three scenarios, and slightly larger for load-independent. The joint angle prediction results of the GP model for all three scenarios were optimal for speed-dependent, with the smallest NRMSE, followed by multi-speed, and worst for speed-independent. The variability of joint motion at different speeds affected the prediction results of the model, and the experimental results also demonstrated that the speed-dependent results were optimal and the speed-independent results were the worst. The average NRMSE of NARX-GP model prediction results for all scenarios was 0.0045 ± 0.0040, which was better than the GP model results (0.1842 ± 0.0602), with a significant difference (p < 0.05).

TABLE 1

Table 1. NRMSE between the estimated joint torque of different models and the measurements (“true” values) of all subjects (mean ± std).

TABLE 2

Table 2. CC between the predicted joint angle of different models and the measurements (“true” values) of all subjects (mean ± std).

3.3.2. Joint Angle Prediction Results for Dataset 2

The proposed method was further tested using the data in dataset 2 to validate the accuracy of the proposed method for predicting the knee joint angle of patients. To reduce the influence of different muscle selections on the prediction results, the sEMG signals of the same pair of muscles, semimembranosus and lateral femoris, were selected for testing in dataset 2 as in dataset 1 and used to build a non-parametric model for joint angle prediction. Dataset 2 included data from the left and right legs of patients with different speeds, so the joint angle prediction results of subjects in three cases of speed-dependent, multi-speed and speed-independent were tested separately for the left and right legs according to the data allocation strategy, and the test results of joint angle prediction for speed of 5 m/s are shown in Figures 11–13, where “NARX-GP” was the joint angle prediction based on the Gaussian process autoregressive model, “GP” was the joint angle prediction based on the Gaussian process model, and “measurement” was the actual measurement of joint angle. The gray shading was the 95% confidence interval for the prediction of the joint angle based on the Gaussian process autoregressive model to describe the uncertainty. From Figures 11–13, it can be concluded that the patient's knee joint NARX-GP model has better angle prediction than the GP model and can achieve good prediction results, and the joint angle prediction can approximate the actual joint angle for both the healthy side and the affected side (right side) of the patient.

FIGURE 11

Figure 11. Joint angle prediction of 5 m/s (speed-dependent).

FIGURE 12

Figure 12. Joint angle prediction of 5 m/s (multi-speed).

FIGURE 13

Figure 13. Joint angle prediction of 5 m/s (speed-independent).

The error assessment and statistical analysis of the prediction results of the NARX-GP model and the GP model for the hemiplegic subject dataset, the NRMSE and CC between the predicted and actual values of joint angles at different speeds are shown in Figures 14, 15. The errors between the predicted results and the actual values of joint angles for both the left and right leg NARX-GP models were small, highly correlated, and significantly better than the GP model.

FIGURE 14

Figure 14. Angle prediction results of NARX-GP model and GP model for hemiplegic subject dataset.

FIGURE 15

Figure 15. Angle prediction results of NARX-GP model and GP model for hemiplegic subject dataset.

Further evaluation and analysis of variance of the prediction results for the left and right legs of the dataset showed that the mean NRMSE of the prediction results of the NARX-GP model for the left and right legs were 0.0063 ± 0.0081 and 0.0032 ± 0.0028, respectively, which were significantly better than those of the GP model (0.1466 ± 0.0127 and 0.1012 ± 0.0092), with significant differences (p < 0.05 and p < 0.05). Joint angle prediction using the NARX-GP model for both the healthy and affected side of the patient was able to have high accuracy with no significant difference (p = 0.1919>0.05). the NRMSE of the NARX-GP model prediction results for all scenarios was 0.0047 ± 0.0063 on average, which was better than the GP model results (0.1239 ± 0.0253) with a significant difference (p < 0.05).

4. Discussion and Conclusion

The EMG signal contains abundant motion information, which is ahead of the actual joint motion, and is often used as a control signal to predict joint motion. Therefore, EMG signals are widely used in applied scientific research related to the development of intelligent rehabilitation technologies and devices. EMG signal based modeling of the neuromusculoskeletal system, as an important component of joint motion prediction, has become a hot topic of research as it is important to help the development of rehabilitation techniques and equipment for patients with sports injuries. There is uncertainty in the neuromusculoskeletal system, and in order for the model to provide a description of the uncertainty, this paper proposes to model the uncertainty using a Gaussian autoregressive model. The muscle activation dynamics model was first introduced into the Gaussian process model to establish a joint angle prediction model based on Gaussian process. Due to the high requirement for joint angle prediction accuracy in practical applications and the fact that the neuromusculoskeletal system is a dynamic non-linear system, the NARX model was introduced into the Gaussian process model to establish a Gaussian autoregressive model to achieve OSA prediction of knee joint angle. The results of different test scenarios on the healthy subjects and hemiplegic subject datasets showed that the designed Gaussian autoregressive model had significantly better prediction accuracy than the Gaussian process model, and there was no significant difference in the prediction accuracy between the affected and healthy sides of the hemiplegic subject, both of which were able to achieve more accurate prediction results for knee angles and could provide uncertainty information.

In this paper, a non-parametric model for knee joint angle prediction was developed from a predictive value-based NARX model approach by mixing a muscle activation kinetic model with a data-driven model. The proposed modeling approach was validated with a publicly available dataset. The proposed method utilizes only the EMG signals of a pair of antagonistic muscles, reducing the cost of EMG signal detection and the complexity of the model. However, there are still some shortcomings in this paper and there are many problems that have not yet been studied with some need for improvement. In this paper, the performance of only one pair of antagonist muscles in the hamstrings and quadriceps was tested, and the effect of sEMG signals from other muscles in the hamstrings and quadriceps as input on the accuracy of knee joint angle prediction can be further tested. In addition, although the knee joint is used as the research object for the study and validation of the model in this paper, the proposed joint angle prediction method is not limited to the knee joint angle prediction. In the subsequent research, the proposed method can be applied to the angle prediction of other joints for relevant testing and validation.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

Conception of the study was conducted by JL, ZS, and WC. FZ processed the data. XC analyzed and interpreted the data with input and support from XC and JL. WC drafted the manuscript. YL, ZS, and WC revised the manuscript critically for important intellectual content. All authors read and approved the manuscript.

Funding

This work was supported in part by the Fujian Province Nature Science Foundation of China under Grant 2019J01544, in part by the National Nature Science Foundation of China under Grant 61773124.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the members of Fuzhou Second Hospital Affiliated to Xiamen University for valuable discussion and thanks to the Fujian Key Laboratory of Medical Instrumentation and Pharmaceutical Technology for providing the experimental environment and experimental equipment.

References

1. Liu J, Kang SH, Xu D, Ren Y, Lee SJ, Zhang LQ. EMG-based continuous and simultaneous estimation of arm kinematics in able-bodied individuals and stroke survivors. Front Neurosci. (2017) 11:480. doi: 10.3389/fnins.2017.00480

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Furukawa Ji, Noda T, Teramae T, Morimoto J. Human movement modeling to detect biosignal sensor failures for myoelectric assistive robot control. IEEE Trans Robot. (2017) 33:846–57. doi: 10.1109/TRO.2017.2683522

CrossRef Full Text | Google Scholar

3. Côté-Allard U, Fall CL, Drouin A, Campeau-Lecours A, Gosselin C, Glette K, et al. Deep learning for electromyographic hand gesture signal classification using transfer learning. IEEE Trans Neural Syst Eng Rehabil. (2019) 27:760–71. doi: 10.1109/TNSRE.2019.2896269

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Corcos DM, Gottlieb GL, Latash ML, Almeida GL, Agarwal GC. Electromechanical delay: an experimental artifact. J Electromyogr Kinesiol. (1992) 2:59–68. doi: 10.1016/1050-6411(92)90017-D

CrossRef Full Text | Google Scholar

5. Faisal AA, Selen LP, Wolpert DM. Noise in the nervous system. Nat Rev Neurosci. (2008) 9:292–303. doi: 10.1038/nrn2258

CrossRef Full Text | Google Scholar

6. Berniker M, Jarc A, Kording K, Tresch M. A probabilistic analysis of muscle force uncertainty for control. IEEE Trans Biomed Eng. (2016) 63:2359–67. doi: 10.1109/TBME.2016.2531083

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Chopp-Hurley JN, Langenderfer JE, Dickerson CR. Probabilistic evaluation of predicted force sensitivity to muscle attachment and glenohumeral stability uncertainty. Ann Biomed Eng. (2014) 42:1867–79. doi: 10.1007/s10439-014-1035-3

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lennart L. System Identification: Theory for the User. Upper Saddle River, NJ: PTR Prentice Hall (1999).

Google Scholar

9. Zeng N, Wang Z, Liu W, Zhang H, Hone K, Liu X. A dynamic neighborhood-based switching particle swarm optimization algorithm. IEEE Trans Cybernet. (2020). doi: 10.1109/TCYB.2020.3029748. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zeng N, Wang Z, Zhang H, Kim KE, Li Y, Liu X. An improved particle filter with a novel hybrid proposal distribution for quantitative analysis of gold immunochromatographic strips. IEEE Trans Nanotechnol. (2019) 18:819–29. doi: 10.1109/TNANO.2019.2932271

CrossRef Full Text | Google Scholar

11. Zeng N, Wang Z, Zineddin B, Li Y, Du M, Xiao L, et al. Image-based quantitative analysis of gold immunochromatographic strip via cellular neural network approach. IEEE Trans Med Imaging. (2014) 33:1129–36. doi: 10.1109/TMI.2014.2305394

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zeng N, Wang Z, Li Y, Du M, Cao J, Liu X. Time series modeling of nano-gold immunochromatographic assay via expectation maximization algorithm. IEEE Trans Biomed Eng. (2013) 60:3418–24. doi: 10.1109/TBME.2013.2260160

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Roberts S, Osborne M, Ebden M, Reece S, Gibson N, Aigrain S. Gaussian processes for time-series modelling. Philos Trans R Soc A Math Phys Sci Eng. (2013) 371:20110550. doi: 10.1098/rsta.2011.0550

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Eslamy M, Alipour K. Synergy-based Gaussian process estimation of ankle angle and torque: conceptualization for high level controlling of active robotic foot prostheses/orthoses. J Biomech Eng. (2019) 141:021002. doi: 10.1115/1.4041767

CrossRef Full Text | Google Scholar

15. Hamaya M, Matsubara T, Noda T, Teramae T, Morimoto J. Learning assistive strategies for exoskeleton robots from user-robot physical interaction. Pattern Recogn Lett. (2017) 99:67–76. doi: 10.1016/j.patrec.2017.04.007

CrossRef Full Text | Google Scholar

16. Kang H, Park F. Motion optimization using Gaussian process dynamical models. Multibody Syst Dyn. (2015) 34:307–25. doi: 10.1007/s11044-014-9441-8

CrossRef Full Text | Google Scholar

17. Schearer EM, Liao YW, Perreault EJ, Tresch MC, Memberg WD, Kirsch RF, et al. Identifying inverse human arm dynamics using a robotic testbed. In: 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. Chicago, IL: IEEE (2014). p. 3585–91. doi: 10.1109/IROS.2014.6943064

CrossRef Full Text | Google Scholar

18. Ullauri JB, Peternel L, Ugurlu B, Yamada Y, Morimoto J. On the EMG-based torque estimation for humans coupled with a force-controlled elbow exoskeleton. In: 2015 International Conference on Advanced Robotics (ICAR). Istanbul: IEEE (2015). p. 302–7. doi: 10.1109/ICAR.2015.7251472

CrossRef Full Text | Google Scholar

19. Xiloyannis M, Gavriel C, Thomik AA, Faisal AA. Gaussian process autoregression for simultaneous proportional multi-modal prosthetic control with natural hand kinematics. IEEE Trans Neural Syst Eng Rehabil. (2017) 25:1785–801. doi: 10.1109/TNSRE.2017.2699598

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yang X, Yan J, Chen Z, Ding H, Liu H. A proportional pattern recognition control scheme for wearable a-mode ultrasound sensing. IEEE Trans Ind Electron. (2019) 67:800–8. doi: 10.1109/TIE.2019.2898614

CrossRef Full Text | Google Scholar

21. Belić JJ, Faisal AA. Decoding of human hand actions to handle missing limbs in neuroprosthetics. Front Comput Neurosci. (2015) 9:27. doi: 10.3389/fncom.2015.00027

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zeng Y, Yang J, Yin Y. Gaussian process-integrated state space model for continuous joint angle prediction from EMG and interactive force in a human-exoskeleton system. Appl Sci. (2019) 9:1711. doi: 10.3390/app9081711

CrossRef Full Text | Google Scholar

23. Raj R, Sivanandan K. Comparative study on estimation of elbow kinematics based on EMG time domain parameters using neural network and ANFIS NARX model. J Intell Syst Fuzzy. (2017) 32:791–805. doi: 10.3233/JIFS-16070

CrossRef Full Text | Google Scholar

24. Li Z, Hayashibe M, Fattal C, Guiraud D. Muscle fatigue tracking with evoked EMG via recurrent neural network: toward personalized neuroprosthetics. IEEE Comput Intell Mag. (2014) 9:38–46. doi: 10.1109/MCI.2014.2307224

CrossRef Full Text | Google Scholar

25. Raj R, Ramakrishna R, Sivanandan KS. A real time surface electromyography signal driven prosthetic hand model using PID controlled DC motor. Biomed Eng Lett. (2016) 6:276–86. doi: 10.1007/s13534-016-0240-4

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gupta R Dhindsa IS Agarwal R. Continuous angular position estimation of human ankle during unconstrained locomotion. Biomed Signal Process Control. (2020) 60:101968. doi: 10.1016/j.bspc.2020.101968

CrossRef Full Text | Google Scholar

27. Liu Z, Wang X, Su M, Le L. Research on rehabilitation training bed with action prediction based on NARX neural network. Int J Imaging Syst Technol. (2019) 29:539–46. doi: 10.1002/ima.22334

CrossRef Full Text | Google Scholar

28. Raj R, Sivanandan K. Elbow joint angle and elbow movement velocity estimation using NARX-multiple layer perceptron neural network model with surface EMG time domain parameters. J Back Rehabil Musculoskel. (2017) 30:515–25. doi: 10.3233/BMR-160525

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zajac FE. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Crit Rev Biomed Eng. (1989) 17:359–411.

PubMed Abstract | Google Scholar

30. Buchanan TS, Lloyd DG, Manal K, Besier TF. Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command. J Appl Biomech. (2004) 20:367–95. doi: 10.1123/jab.20.4.367

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Lloyd DG, Besier TF. An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo. J Biomech. (2003) 36:765–76. doi: 10.1016/S0021-9290(03)00010-1

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Li Y, Chen W, Yang H, Li J, Zheng N. Joint torque closed-loop estimation using NARX neural network based on sEMG signals. IEEE Access. (2020) 8:213636–46. doi: 10.1109/ACCESS.2020.3039983

CrossRef Full Text | Google Scholar

33. Schulz E, Speekenbrink M, Krause A. A tutorial on Gaussian process regression: modelling, exploring, and exploiting functions. J Math Psychol. (2018) 85:1–16. doi: 10.1016/j.jmp.2018.03.001

CrossRef Full Text | Google Scholar

34. Williams CKI. Gaussian Processes for Machine Learning. Taylor & Francis Group (2006).

Google Scholar

35. Stein ML. Interpolation of Spatial Data: Some Theory for Kriging. Springer Science & Business Media (2012).

Google Scholar

36. Leontaritis I, Billings SA. Input-output parametric models for non-linear systems part I: deterministic non-linear systems. Int J Control. (1985) 41:303–28. doi: 10.1080/0020718508961129

CrossRef Full Text | Google Scholar

37. Leontaritis I, Billings SA. Input-output parametric models for non-linear systems part II: stochastic non-linear systems. Int J Control. (1985) 41:329–44. doi: 10.1080/0020718508961130

CrossRef Full Text | Google Scholar

38. Hamner SR, Delp SL. Muscle contributions to fore-aft and vertical body mass center accelerations over a range of running speeds. J Biomech. (2013) 46:780–7. doi: 10.1016/j.jbiomech.2012.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Meyer AJ, Patten C, Fregly BJ. Lower extremity EMG-driven modeling of walking with automated adjustment of musculoskeletal geometry. PLoS ONE. (2017) 12:e0179698. doi: 10.1371/journal.pone.0179698

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Babak S. Biostatistics With R: An Introduction to Statistics Through Biological Data. New York, NY: Springer New York (2012).

Google Scholar