US12106204B2 - Adaptive brain-computer interface decoding method based on multi-model dynamic integration - Google Patents

Abstract

The present invention discloses an adaptive brain-computer interface decoding method based on multi-model dynamic ensemble, where a traditional state-space model is improved, and a set of measurement functions instead of one fixed measurement function are used to dynamically characterize a relationship between observation variables and state variables; and, by using a pool of linear and nonlinear decoders, and in a decoding process of a brain-computer interface system, decoders are automatically switched according to the data, so as to realize adaptive brain signal decoding. Through the above multi-model ensemble strategy, linear and nonlinear decoder capabilities can be integrated, the accuracy and stability of the brain-computer interface system can be improved, and decoding unstability caused by the non-stationary neural signal of the brain-computer interface system can be solved to a certain extent.

Description

This is a U.S. national stage application of PCT Application No. PCT/CN2021/126580 under 35 U.S.C. 371, filed Oct. 27, 2021 in Chinese, claiming priority to Chinese Patent Applications No. 202011593629.3, filed Dec. 29, 2020, all of which are hereby incorporated by reference.

FIELD OF TECHNOLOGY

The present invention relates to the field of motor neural signal decoding, in particular to an adaptive brain-computer interface decoding method based on multi-model dynamic ensemble.

BACKGROUND TECHNOLOGY

Brain-computer interfaces directly decode motion intentions from neural signals for external device control. Intracortical brain-computer interfaces use neural signals recorded by implanted electrode arrays to extract motion-related information. Advances in the intracortical brain-computer interfaces have led to the significant development of exoskeleton devices and computer cursor control.

In the intracortical brain-computer interface systems, neural decoding algorithms play an important role. Scholars have proposed many algorithms to decode motion information from the neural signals, including population vector methods, optimal linear estimation methods, deep neural networks, and recursive Bayesian decoders. Among these methods, the Kalman filter combines the evolution process of a trajectory as prior knowledge to make prediction more accurate, so it is widely used in online cursor control and exoskeleton control, achieving a best performance so far.

For example, Chinese patent document with publication number CN107669416A discloses a wheelchair system and control method based on continuous-brisk motor imagery neural decoding, comprising an electroencephalogram signal acquisition system, a decoding device and an electric wheelchair connected in sequence; the decoding device comprises a continuous motor imagery electroencephalogram signal processing module and a brisk motor imagery electroencephalogram signal processing module that are arranged in series, where the two modules are respectively used to decode the continuous motor imagery electroencephalogram signal and the brisk motor imagery electroencephalogram signal transmitted by the electroencephalogram signal acquisition system.

However, the unstability of the neural signals is a great challenge in the neural decoding process. Most of the neural decoders used in current intracortical brain-computer interfaces assume a stable functional relationship between the neural signals and kinematics, and therefore use a fixed decoding model. However, during online decoding, signals generated by neurons occasionally introduce noise or disappear altogether; at the same time, because of the plasticity of the neurons, brain activity may also change over time. Because of the noise and change described above, a neural signal-to-kinematic mapping function may be unstable and change continuously over time. A fixed decoding function can lead to an unstable and inaccurate decoding result, so it needs to be retrained every time period to maintain a certain performance.

Decoders proposed for coping with neural variabilities can be classified into two categories. The first category of the decoders, which are the most common decoders today, still uses fixed models, relying on periodic retraining or incrementally updating model parameters online to maintain performance. The second category of the decoders uses dynamic models to track the changes in the neural signals, which can avoid the cost of retraining and are more suitable for long-term decoding tasks. However, few studies in such methods dare to directly model unstable neural signals.

Therefore, in motor neural decoding, how to use dynamic models to model the unstability of the neural signals to achieve a stable decoding performance is an important problem to be solved in the current motor neural decoding field.

SUMMARY OF THE INVENTION

The present invention provides an adaptive brain-computer interface decoding method based on multi-model dynamic ensemble, which can greatly reduce the influence of unstability of the neural signals on decoding performance and improve the robustness of a decoder.

An adaptive brain-computer interface decoding method based on multi-model dynamic ensemble, where the method comprises the following steps:

(1) obtaining the original motor neural signal from electrodes, and after pre-processing the original motor neural signal, dividing the signal into a training set, a validation set and a test set;

(2) building a state-space decoder based on the multi-model dynamic ensemble, which is specified as follows:

• • (2-1) with an improved state-space model, using a pool of candidate models to dynamically characterize a relationship between observation variables and state variables, where the pool of candidate models comprises a linear function, a second-order polynomial function, and two neural networks. Here, the observation variables are neural signals, and the state variables are motion signals;
  • (2-2) dynamically combining the candidate models according to a Bayesian update mechanism as an observation function of the state-space model;
  • (3) using the training set and the validation set to train and evaluate the improved state-space model, and obtaining the parameters of different candidate models by training;
  • (4) using the test set to test the performance and robustness of the model.

In the method of the present invention, the observation function in a Kalman filter is improved and replaced with a pool of models that can be dynamically assembled online, so as to better adapt to the changes of the neural signals. During the decoding process, these models are automatically selected and combined according to the Bayesian update mechanism, thereby greatly reducing the influence of the unstability of the neural signals on the decoding performance and improving the robustness of the decoder.

In the step (1), the specific process of the pre-processing is as follows:

• • selecting an appropriate window size to calculate the firing rate of the neural signal, intercepting a valid data segment according to a state label, normalizing and smoothing the data respectively, and obtaining the pre-processed motor neural signal.

Preferably, the data in paradigm preparation and return phases can be removed, and an actual operation phase can be selected for analysis.

Preferably, when performing the normalizing and smoothing operations, a z-score function and a filter function can be used to normalize and smooth each dimension of the motion signal, and smooth the firing rate of each neuron of the neural signal. The specific smoothing coefficients can be selected according to actual needs.

In the step (2), the expression of the improved state-space model is as follows:
x _k=ƒ(x _k-1)+v _k-1
y _k =h _{H k} (x _k)+n _k
where, k represents a discrete time step; x_k∈R^{d x} represents a state of interest; ƒ(·) represents a state transition function; y_k ∈R^{d y} represents an observation or measurement variable; n_k represents an independent and identically distributed observation noise;

∈{1, 2, . . . , M} represents a model index in the observation function h(·),

=m represents that the model in action at the moment k is h_m.

In the step (2), the Bayesian update mechanism is as follows:

dynamically selecting and switching among a pool of candidate models in a data-driven manner, where given an observation sequence y_0:k, the posterior probability of the state at the moment k is:
p(x _k |y _0:k)=Σ_m=1 ^M p(x _k |

=m,y _0:k)p(

=m|y _0:k)

Where, p(x_k|

=m, y_0:k) is the posterior of the state when

=m; p(

=m|y_0:k) is the posterior probability of selecting the m-th model; M represents the number of the candidate models.

Considering that p(

=m|y_0:k) is obtained by iteration of p(

=m|y_0:k-1), preferably, the Bayesian formula can be used to calculate the posterior probability of selecting the m-th model at the moment k. The formula thereof is as follows:

$p\left(=m|y_{0:k}\right)=\frac{p\left(=m|y_{0:k-1}\right)p_m\left(y_k|y_{0:k-1}\right)}{{\Sigma }_{j=1}^{M}p\left(=j|y_{0:k-1}\right)p_j\left(y_k|y_{0:k-1}\right)}$

Where, p(

=M|y_0:k-1) is a priorprior probability of selecting the m-th model; p_m(y_k|y_0:k-1) is a marginal likelihood of selecting the m-th model.

Preferably, a forgetting factor α can be determined to retain part of the historical information, that is, the prior probability of selecting the m-th model at the current moment is influenced by the posterior probability of selecting the m-th model at the previous moment. The priorprior probability of selecting the m-th model is calculated as follows:

$p\left(=m|y_{0:k-1}\right)=\frac{{p\left(=m|y_{0:k-1}\right)}^{\alpha }}{{\Sigma }_{j=1}^M{p\left(=j|y_{0:k-1}\right)}^{\alpha }}$

where, p(

=m|y_0:k-1) is the probability of selecting the m-th model at the k-1-th moment; a (0<α<1) is a forgetting factor.

Preferably, at the moment k, the calculation method of the marginal likelihood of selecting the m-th model is as follows:
p _m(y _k |y _0:k-1)=∫p _m(y _k |x _k)p(x _k |y _0:k-1)dx _k

where p_m(y_k|x_k) is a likelihood function for the m-th model.

In step (4), adopting a particle filter algorithm for the state estimation, and calculating p(x_k|

=m, y_0:k) and p(

=m|y_0:k) based on particles.

Based on the traditional state-space model and particle filtering, the present invention proposes a multi-model dynamic ensemble method applied to the neural signal decoding, which reduces the influence caused by the unstability of the neural signal to a certain extent. It outperforms the well-known Kalman filter algorithm in the test, proving the effectiveness of the method.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic figure of a multi-model dynamic ensemble process in the method of the present invention;

FIG. 2 is a comparison chart of the success rate of a Kalman filter and multi-model dynamic ensemble for five consecutive days in an online cursor control experiment in an embodiment of the present invention;

FIG. 3 is a trajectory comparison figure of a Kalman filter and multi-model dynamic ensemble in an online cursor control experiment on Dec. 14, 2020, according to an embodiment of the present invention.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be pointed out that the following embodiments are intended to facilitate the understanding of the present invention, but do not have any limiting effect on it.

All clinical and experimental procedures in this embodiment were approved by the Medical Ethics Committee of the Second Affiliated Hospital of Zhejiang University (ethical review number 2019-158). Informed consent was obtained orally by the participants and their immediate family members and signed by their legal representatives. The volunteer was a 74-year-old male who suffered a car accident three years ago and was paralyzed after a C4 cervical spine trauma. The volunteer is only able to move the parts above the neck and has normal language skills and comprehension ability of all tasks. For a limb motion behavior, the patient has a skeletal muscle strength score of 0/5 and a complete loss of limb motion control.

Two 96-channel Utah intracortical microelectrode arrays (4 mm×4 mm, 1.4 mm long Utah arrays, Blackrock Microsystems, Salt Lake City, Utah, USA) were implanted in the left primary motor cortex of the volunteer, and the other was located in the central hand knob area 2 mm away. Electron Computed Tomography and functional MRI were used to guide the implantation. During the implantation, the participant was asked to perform flexion and extension motions of the elbow, and the scanning of the functional MRI was used to confirm an activation area in the motor cortex. The operation used a robotic arm to assist in implanting electrodes during the surgery. The participant had a week of rehabilitation before starting the neural signal recording and BCI training tasks.

The participant performed the BCI training tasks each weekday and rested on weekends. The training time was approximately 3 hours per day and comprised of preparation for signal recording, impedance testing, spike classification and paradigm tasks. The experiment was stopped once the participant felt tired or developed physical abnormalities, such as a fever or a urinary tract infection. The entire experiment lasted half a year, recording the neural signals for 12 weeks.

A 2D cursor control task was performed on a computer monitor 1.5 m in front of the volunteer. The task contained eight targets at the top, bottom, left, right and four corners of the screen. The targets can be displayed sequentially in a clockwise or pseudo-random fashion. The relevant parameters of the task can be configured in task settings, comprising a distance from the target to the center of the screen, the diameter of a target ball, the maximum trial time, an arrival threshold, the minimum dwell time and the number of repetitions after failure. For each experimental trial, the volunteer was asked to control a blue cursor ball from the center of the screen to a red target ball. In a successful trial, the distance between the centers of the blue and red balls should be less than a preset arrival threshold, and dwell time should not be shorter than a preset dwell time. The maximum duration of each trial is set from 2 to 5 seconds. An audible cue is given at the end of the trial to indicate its success. A maximum of 4 attempts were made after each trial failed. The default setting was to display the target balls at 15 cm from the center of the screen, with a diameter of 5 cm.

For the adaptive brain-computer interface decoding method based on multi-model dynamic ensemble proposed by the present invention, the specific implementation was as follows:

Step 1, motor neural signal pre-processing: obtaining an original motor neural signal recorded from the electrodes, selecting an appropriate window size to calculate the firing rate of the neural signal, intercepting a valid data segment according to a state label, normalizing and smoothing the data, and obtaining the pre-processed motor neural signal; dividing the data into a training set, a validation set, and a test set according to a proper proportion;

Specifically, the experiment recorded the neural signals by using a Neuroport system (NSP, Blackrock Microsystems). The neural activity is amplified, digitized, and recorded at 30 KHz. Thresholds for neural action potential detection were set from −6.5 RMS to −5.5 RMS for each high-pass filtered (250 Hz cutoff frequency) electrode by using a Central software package (Blackrock Microsystem). At the beginning of each routine task, the researchers manually sort the spikes, which takes about 25 to 35 minutes. A peak activity was converted to firing rate in 20 ms and low-pass filtered by using an exponential smoothing function with a 450 ms window.

Step 2, a state-space model based on the dynamic ensemble:

2-1. An improved observation model: varying a traditional state-space model, using a set of functions (i.e., candidate models) instead of a fixed function to dynamically characterize the relationship between observation variables and state variables;

The traditional state-space model consists of a state transition function ƒ(·) and an observation function h(·):
x _k=ƒ(x _k-1)+v _k-1
y _k =h(x _k)+n _k

Where, k represents a discrete time step; x_k∈R^{d x} represents the state that is concerned about; y_k∈R^{d y} represents the observation variable; v_k, n_k represents a state transition noise and an observation noise which are independent and identically distributed.

In the field of the neural decoding, states and observations represent motion trajectories and neural signals, respectively. Given a neural signal sequence y_0:k, the task of the neural decoding is to iteratively estimate the probability density of x_k. The Kalman filter can provide an analytically optimal solution when both the neural signal and the motion signal conform to a linear Gaussian assumption.

The observation function h(·) of the traditional state-space model is determined in advance and cannot adapt to changing neural signals. The observation model in step 2 improves upon it by allowing the observation function to be adjusted online. The improved observation function is expressed as follows:
y _k=

(x _k)+n _k

Where, k represents a discrete time step; x_k∈ R^{d x} represents a state of interest; y_k∈R^{d y} represents the observation variable; n_k represents an independent and identically distributed observation noise;

∈{1, 2, . . . , M} represents a model index in the observation equation. In particular,

=m represents the model in action at the moment k is h_m.

The models in a model set

, similar to neural encoders, can convert input kinematic parameters x_t into the neural signals y_t. In the multi-model dynamic ensemble method, we learn a pool of forms of encoding models, such as linear functions, polynomial functions, and neural networks, to enhance the ensemble effectiveness of online models. Past research has demonstrated an excellent decoding ability of linear models, while non-linear models have strong robustness in the presence of noise. As shown in FIG. 1 , the multi-model dynamic ensemble method demonstrates the advantages of different models to deal with unstable neural signals.

Specifically, the model pool of this example contains four models. The first model is a linear function H_l(x)=Wx+b, which is similar to the observation model in the Kalman filter. The second model is a second-order polynomial function H_p (x)=a₀+a₁x+a₂x². The third model H_nn1(·) and the fourth model H_nn2 (·) are two neural networks with different parameter sizes. Both the neural networks contain one hidden layer. The number of neurons in the input layer is equal to the dimensions of the motion parameters x_t, and the number of neurons in the output layer is equal to the dimensions of the neural signal y_t. The number of neurons in the hidden layers are 30 and 50, respectively.

For H_l(·) and H_p (·), the parameters are estimated by using the least square method. For H_nn1(·) and H_nn2 (·), the parameters are optimized by using the Adam algorithm, a learning rate is set to 0.01, and a weight decay is set to 1e-4. At the same time, an early-stopping method is adopted to alleviate the overfitting problem.

2-2 Dynamically combining the candidate models according to a Bayesian update mechanism as an observation function of the state-space model;

In step 2, a set of candidate models are used to characterize the functional relationship between the observation signal and the state to be predicted. The Bayesian update mechanism dynamically selects dominant models among these models in a data-driven manner. Given an observation sequence y_0:k, the posterior probability of the state at the moment k is:
p(x _k |y _0:k)=Σ_m=1 ^M p(x _k |

=m,y _0:k)p(

=m|y _0:k)

Where, p(x_k|H_k=m, y_0:k) is the posterior of the state when

=p(

=m|y_0:k) is the posterior probability of selecting the m-th model.

Considering that p(

=m|y_0:k) is obtained by iteration of p(

=m|y_0:k-1), preferably, the Bayesian formula can be used to calculate the posterior probability of selecting the m-th model at the moment k.

$p\left(=m|y_{0:k}\right)=\frac{p\left(=m|y_{0:k-1}\right)p_m\left(y_k|y_{0:k-1}\right)}{{\Sigma }_{j=1}^{M}p\left(=j|y_{0:k-1}\right)p_j\left(y_k|y_{0:k-1}\right)}$

Where, p(

=m|y_0:k-1) is a priorprior probability of selecting the m-th model; p_m(y_k|y_0:k-1) is a marginal likelihood of selecting the m-th model.

Preferably, a forgetting factor α can be determined to retain part of the historical information, that is, the prior probability of selecting the m-th model at the current moment is influenced by the posterior probability of selecting the m-th model at the previous moment:

$p\left(=m|y_{0:k-1}\right)=\frac{{p\left(=m|y_{0:k-1}\right)}^{\alpha }}{{\Sigma }_{j=1}^M{p\left(=j|y_{0:k-1}\right)}^{\alpha }}$

Where, p(

=m|y_0:k-1) is the probability of selecting the m-th model at the k-1-th moment; α(0<α<1) is a forgetting factor.

In this example, the forgetting factor is set from 0.1 to 0.9. When the forgetting factor is set to 0.9, the influence of the previous kinematic state is considered to a great extent, so the decoded kinematic state is smoother.

Preferably, at the moment k, the calculation method of the marginal likelihood of selecting the m-th model is as follows:
p _m(y _k |y _0:k-1)=∫p _m(y _k |x _k)P(x _k |y _0:k-1)dx _k

Where p_m (y_k|x_k) is a likelihood function for the m-th model.

Step 3, candidate model learning and state estimation: using the training set and the validation set to obtain a set of different candidate models;

Preferably, based on the particle filter algorithm, estimating the kinematic state for each candidate model by using a weighted particle set at each time step, and calculating p(x_k|

=m, y_0:k) and p(

=m|y_0:k) in step 2.

Assuming that this moment is the beginning of the k-th time step, and p

=m|y_0:k), m=1, . . . , M and the weight of the particle set are known, where N_S represents the size of the particle set and x_k-1 ⁱ is the particle with the weight of ω_k-1 ⁱ. Assuming p(x_k-1|y_0:k-1)≅Σ_i=1 ^{N s} ω_k-1 ⁱδ(x_k-1−x_k-1 ⁱ), where δ(·) represents a Dirac δ function, here is how to compute p(x_k|

=m, y_0:k) and p(

=m|y_0:k) by using particle filtering:

3-1, p(x_k|

=m, y_0:k) based on the particle filtering. First, x_k ⁱ, is obtained from the state transition prior p(x_k|x_k-1 ⁱ),i=1, . . . , N_S, and then according to the principle of importance sampling, we can get:
p(x _k |

=y _0:k)≈Σ_i=1 ^{N s} ω_m,k ⁱδ(x _k −x _k ⁱ)

Where, ω_m,k ⁱ∝ω_k-1 ⁱp_m(y_k|x_k ⁱ), Σ_i=1 ^{N s} ω_m,k ⁱ=1. ω_m,k ⁱ represents the normalized importance weight of the i-th particle when the m-th model is selected.

In this example, the training set is used as the particle set to avoid the unstability caused by randomly scattering the particles.

3-2, p(

=m|y_0:k) based on the particle filtering. Given p(

=m|y_0:k-1), first, the priorprior probability of

=m is computed by using the forgetting factor in step 2. The posterior probability p_m(y_k=m|y_0:k) of

=m can be calculated when the likelihood p_m(y_k|y_0:k-1), m=1, . . . , M is known. In the step 3-1, the state transition prior is used as the importance function, i.e., q(x_k|x_k-1, y_0:k)=p(x_k|x_k-1). Therefore, the distribution of x_k can be approximated by the particles, i.e, p(x_k|y_0:k-1)≈Σ_i=1 ⁱδ_{x k} ⁱ. Further, we get:
p _m(y _k |y _0:k-1)≈Σ_i=1 ^{N s} ω_k-1 ⁱ p _m(y _k |x _k ⁱ)

Note that the particle filtering usually suffers from particle degradation, and after a few iterations, only a small number of particles have high weights. Therefore, we adopt a resampling method to remove particles with too small or too large weights to alleviate the particle degradation problem.

Step 4, performance evaluation of the multi-model dynamic ensemble method: using this method to compare with other methods in the test set data to evaluate the performance and robustness of this method.

This example was an online experiment with closed-loop calibration. Two phases of calibration and testing were contained in each experiment. Each session contains several blocks, and each block consists of 16 trials, where 2 trials for each target. The neural decoder was trained during the calibration phase. There were two observation blocks in the calibration phase, and the computer automatically completed the cursor control task. The neural signals of the volunteer during the observation phase were used to initialize the decoder. Subsequently, the decoder was further adjusted by using computer-assistant and ortho-impedance assistance, and the assistant ratio decreased with the increase of the number of blocks. The calibration process lasted about 10 minutes, and then the neural signals coming in the test phase were decoded by using the trained decoder.

In order to compare and illustrate the effectiveness of the neural signal decoding method proposed by the present invention, the success rate statistics were performed on the online experimental data for five consecutive days, and FIG. 2 was obtained. Compared with the Kalman filter, the online success rate of the multi-model dynamic ensemble method was increased by 2.30% (the big ball), 5.68% (the small ball), 13.85% (tracing), and the arrival time was shortened by 3.52% (the big ball), 8.81% (the small ball), 13.46% (tracing), and the performance of different days was more stable, the variance was 47.50% smaller than the Kalman filter (the big ball), 70.42% (the small ball), 64.51% (tracing). FIG. 3 is the trajectory diagram of the ball test phase on Dec. 14, 2020. It can be seen that the trajectory solved by the Kalman filter has a strong bias and cannot reach each target point on time.

However, the trajectory distribution solved by the multi-model dynamic ensemble method is relatively uniform and successfully reaches each target point. It can be seen from the above results that the performance of the multi-model dynamic ensemble method is more stable.

The above-mentioned embodiments describe the technical solutions and beneficial effects of the present invention in detail. It should be understood that the above-mentioned embodiments are only specific embodiments of the present invention, and are not intended to limit the present invention. Any modifications, additions and equivalent replacements made within the scope of the principles of the present invention shall be included within the protection scope of the present invention.

Claims (8)

The invention claimed is:

1. A computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor, where the method comprises the following steps:

(1) obtaining an original motor and neural signal collected by electrodes implanted in a brain of a subject, specifically comprising: planting the electrodes in a primary motor cortex of the subject for recording the original motor and neural signal; the subject being asked to imagine performing a two dimensional (2D) cursor control task on a computer monitor to move a cursor from a center of the computer monitor through a trajectory to a target; recording the subject's neural signals simultaneously during the cursor control task to obtain the original motor and neural signal; and after pre-processing the original motor and neural signal, dividing the signal into a training set, a validation set and a test set;

(2) building a state-space decoder based on the multi-model dynamic ensemble, which is specified as follows:

(2-1) with a state-space model, using a pool of candidate models to dynamically characterize a relationship between observation variables and state variables, where the pool of candidate models comprise a linear function H₁(x)=Wx+b, wherein H_l(x) represents a linear mapping function from movement states x to neural signals, W and b represent linear regression parameters; a second-order polynomial function H_p(x)=a₀+a₁x+a₂x², wherein H_p(x) represents a polynomial mapping function from movement states x to neural signals, where a₀, a₁ and a₂ represent the polynomial regression parameters for the term of bias, x and x², respectively, and two neural networks H_nn1(·) and H_nn2(·), and, the observation variables are neural signals, and the state variables are kinematic signals;

(2-2) dynamically combining the candidate models according to a Bayesian update mechanism as an observation function of the state-space model;

(3) using the training set and the validation set to train and evaluate the state-space model, and obtaining parameters of different candidate models by training;

(4) using the test set to test a performance and robustness of the state-space model, applying a neural signal decoding, and obtaining a state estimation after inputting the pre-trained neural signal to be decoded;

where in the step (2-1), an expression of the state-space model is as follows:

X _k =f(x _k-1)+V _k-1

y _k =h _{H k} (x _k)+n _k

where, k represents a discrete time step; x_k∈R^{d x} represents a state of interest; ƒ(·) represents a state transition function; y_k ∈R^{d y} represents an observation or measurement variable; n_k represents an independent and identically distributed observation noise; H_k ∈{1, 2, . . . , M} represents a model index in the observation function h(·), H_k=m represents that the model used at the moment k is h_m;

for H_l(·) and H_p(·), the parameters are estimated by using the least square method, for H_nn1(·) and H_nn2(·), the parameters are optimized by using an Adam algorithm, a learning rate is set to 0.01, and a weight decay is set to 1e-4; at a same time, an early-stopping method is adopted to alleviate an overfitting problem;

where in the step (2-2), the Bayesian update mechanism is as follows:

dynamically selecting and switching among a pool of candidate models in a data-driven manner, where given an observation sequence y_0:k, a posterior probability of a state at the moment k is:

p(x _k |y _0:k)=Σ_m=1 ^M p(x _k |

=m,y _0:k)p(

=m|y _0:k)

where, p(x_k|H_k=m, y_0:k) is the posterior of the state when H_k=m; p(H_k=m|y_0:k) is the posterior probability of selecting the m-th model; M represents the number of the candidate models;

wherein, the step (3) specifically comprises the following steps:

step (3.1): based on a particle filter algorithm, estimating motion state for each candidate model by using a weighted particle set at each time step, and calculating p(x_k|H_k=m, y_0:k) and p(H_k=m|y_0:k) in the step (2);

step (3.2): assuming that this moment is beginning of the k-th time step, and p(H_k-1=m|y_0:k), m=1, . . . , M and a particle set with a weight are known, wherein N_s represents the size of the particle set and x_k-1 ⁱ is the particle with the weight of ω_k-1 ⁱ, and assuming p(x_k-1|y_0:k-1)≅Σ_i=1 ^{N s} ω_k-1 ⁱδ(x_k-1−x_k-1 ⁱ), wherein δ(·) represents a Dirac δ function, calculating p(x_k|

=m, y_0:k) and p(

=m|y_0:k) by using particle filtering as follows:

step (3.2.1): obtaining x_k ⁱ from the state transition priori p(x_k|X_k-1 ⁱ), i=1, . . . , N_s, and then according to a principle of importance sampling, obtaining:

p(x _k |H _k=m,y_0:k)≈Σ_i=1 ^{N s} ω_m,k ⁱδ(x _k −x _k ⁱ)

wherein, ω_m,k ⁱ∝ω_k-1 ⁱp_m(y_k|x_k ⁱ), Σ_i=1 ^{N s} ω_m,k ⁱ=1·ω_m,k ⁱ represents the normalized importance weight of the i-th particle when the m-th model is selected;

step (3.2.2): assuming p(H_k-1=m|y_0:k-1), first, computing the priori probability of H_k=m by using the forgetting factor in the step (2), calculating posterior probability p(H_k=m|y_0:k) of H_k=m when the likelihood P_m(y_k|y_0:k-1), m=1, . . . , M is known; in the step (3.2.1), using state transition priori is the function, q(x_k|x_k-1, y_0:k)=p(x_k|x_k-1), approximating the distribution of x_k by the particles, p(x_k-1|y_0:k-1)≅Σ_i=1 ^{N s} ω_k-1 ⁱδx^k _i, and obtaining:

p _m(y _k |y _0:k-1)≈Σ_i=1 ^{N s} ω_k-1 ⁱ p _m(y _k |x _k ⁱ)

step (5): outputting decoded neural signal to an exoskeleton device or a computer cursor control device.

2. The computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor according to claim 1, where in step (1), the specific process of the pre-processing is as follows:

selecting an appropriate window size to calculate a firing rate of the neural signal, intercepting a valid data segment according to a state label, normalizing and smoothing the data respectively, and obtaining the pre-processed motor and neural signal.

3. The computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor according to claim 1, where the posterior probability of selecting the m-th model is calculated as follows:

$p\left(=m|y_{0:k}\right)=\frac{p\left(=m|y_{0:k-1}\right)p_m\left(y_k|y_{0:k-1}\right)}{{\Sigma }_{j=1}^{M}p\left(=j|y_{0:k-1}\right)p_j\left(y_k|y_{0:k-1}\right)}$

where, p(H_k=m|y_0:k-1) is a prior probability of selecting the m-th model; P_m (y_k|y_0:k-1) is a marginal likelihood of selecting the m-th model.

4. The computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor according to claim 3, where the prior probability of selecting the m-th model is calculated as follows:

$p\left(=m|y_{0:k-1}\right)=\frac{{p\left(=m|y_{0:k-1}\right)}^{\alpha }}{{\Sigma }_{j=1}^M{p\left(=j|y_{0:k-1}\right)}^{\alpha }}$

where, p(H_k-1=m|y_0:k-1) is the probability of selecting the m-th model at the k-1-th moment; α is a forgetting factor, 0<α<1.

5. The computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor according to claim 3, where the marginal likelihood of selecting the m-th model is calculated as follows:

p _m(y _k |y _0:k-1)=∫p _m(y _k |x _k)P(x _k |y _0:k-1)dx _k

where p_m(y_k|x_k) is a likelihood function for the m-th model.

6. The computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor according to claim 3, where in the step (4), adopting a particle filter algorithm for the state estimation, and calculating p(x_k|H_k=m,y_0:k) and p(H_k=m|y_0:k) based on particles.

7. A computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor, where the method comprises the following steps:

(1) obtaining an original motor and neural signal collected by electrodes implanted in a brain of a subject, specifically comprising: planting the electrodes in a primary motor cortex of the subject for recording the original motor and neural signal; the subject being asked to imagine performing a two dimensional (2D) cursor control task on a computer monitor to move a cursor from a center of the computer monitor through a trajectory to a target; recording the subject's neural signals simultaneously during the cursor control task to obtain the original motor and neural signal; and after pre-processing the original motor and neural signal, dividing the signal into a training set, a validation set and a test set;

(2) building a state-space decoder based on the multi-model dynamic ensemble, which is specified as follows:

(2-1) with a state-space model, using a pool of candidate models to dynamically characterize a relationship between observation variables and state variables, where the pool of candidate models comprise a linear function, a second-order polynomial function, and two neural networks, and, the observation variables are neural signals, and the state variables are kinematic signals;

(2-2) dynamically combining the candidate models according to a Bayesian update mechanism as an observation function of the state-space model;

(3) using the training set and the validation set to train and evaluate the state-space model, and obtaining parameters of different candidate models by training;

(4) using the test set to test a performance and robustness of the state-space model, applying a neural signal decoding, and obtaining a state estimation after inputting the pre-trained neural signal to be decoded;

step (5): outputting decoded neural signal to an exoskeleton device or a computer cursor control device;

where in step (1), the specific process of the pre-processing is as follows:

selecting an appropriate window size to calculate a firing rate of the neural signal, intercepting a valid data segment according to a state label, normalizing and smoothing the data respectively, and obtaining the pre-processed motor and neural signal;

where in the step (2), an expression of the state-space model is as follows:

X _k =f(x _k-1)+V _k-1

y _k =h _{H k} (x _k)+n _k

where, k represents a discrete time step; x_k∈ R^{d x} represents a state of interest; f(·) represents a state transition function; y_k∈R^{d y} represents a state of observation or measurement variable; n_k represents an independent and identically distributed observation noise; H_k ∈{1,2, . . . , M} represents a model index in the observation function h(·), H_k=m represents that the model used at the moment k is h_m;

where in the step (2), the Bayesian update mechanism is as follows:

dynamically selecting and switching among a pool of candidate models in a data-driven manner, where given an observation sequence y_0:k, a posterior probability of a state at the moment k is:

p(x _k |y _0:k)=Σ_m=1 ^M p(x _k |H _k =m,y _0:k)p(H _k =m|y _0:k)

where, p(x_k|H_k=m,y_0:k) is the posterior of the state when H_k=m; p(H_k=m|y_0:k)is the posterior probability of selecting the m-th model; M represents the number of the candidate models;

where the posterior probability of selecting the m-th model is calculated as follows:

$p\left({\mathscr{H}}_k=m❘y_{0:k}\right)=\frac{p\left({\mathscr{H}}_k=m❘y_{0:k-1}\right)p_m\left(y_k❘y_{0:k-1}\right)}{{\sum }_{j=1}^{M}p\left({\mathscr{H}}_k=j❘y_{0:k-1}\right)p_j\left(y_k❘y_{0:k-1}\right)}$

where, p(H_k=m|y_0:k-1) is a prior probability of selecting the m-th model; P_m (y_k|y_0:k-1) is a marginal likelihood of selecting the m-th model;

where the prior probability of selecting the m-th model is calculated as follows:

$p\left({\mathscr{H}}_k=m❘y_{0:k-1}\right)=\frac{{p\left({\mathscr{H}}_{k-1}=m❘y_{0:k-1}\right)}^{\alpha }}{{\sum }_{j=1}^M{p\left({\mathscr{H}}_{k-1}=j❘y_{0:k-1}\right)}^{\alpha }}$

where, p(H_k-1=m|y_0:k-1) is the probability of selecting the m-th model at the k-1-th moment; α is a forgetting factor, 0<α<1;

where the marginal likelihood of selecting the m-th model is calculated as follows:

p _m(y _k |y _0:k-1)=∫p _m(y _k |x _k)P(x _k |y _0:k-1)dx _k

where p_m (y_k|x_k) is a likelihood function for the m-th model;

wherein, the step (3) specifically comprises the following steps:

step (3.1): based on a particle filter algorithm, estimating motion state for each candidate model by using a weighted particle set at each time step, and calculating p(x_k|H_k=m, y_0:k) and p(H_k=m|y_0:k) in the step (2);

step (3.2): assuming that this moment is beginning of the k-th time step, and p(H_k-1=m|y_0:k), m=1, . . . , M and a particle set with a weight are known, wherein N_s represents the size of the particle set and x_k-1 ⁱ is the particle with the weight of ω_k-1 ⁱ, and assuming p(x_k-1|y_0:k-1)≅Σ_i=1 ^{N s} ω_k-1 ⁱδ(x_k-1−x_k-1 ⁱ), wherein δ(·) represents a Dirac δ function, calculating p(x_k|

=m, y_0:k) and p(H_k=m|y_0:k) by using particle filtering as follows:

step (3.2.1): obtaining x_k ⁱ from the state transition priori p(x_k|X_k-1 ⁱ), i=1, . . . , N_s, and then according to the principle of importance sampling, obtaining:

p(x _k |H _k=m,y_0:k)≈Σ_i=1 ^{N s} ω_m,k ⁱδ(x _k −x _k ⁱ)

wherein, ω_m,k ⁱ∝ω_k-1 ⁱp_m(y_k|x_k ⁱ), Σ_i=1 ^{N s} ω_m,k ⁱ=1·ω_m,k ⁱ represents the normalized importance weight of the i-th particle when the m-th model is selected;

step (3.2.2): assuming p(H_k-1=m|y_0:k-1), first, computing the priori probability of H_k=m by using the forgetting factor in the step (2), calculating posterior probability p(H_k=m|y_0:k) of H_k=m when the likelihood P_m(y_k|y_0:k-1), m=1, . . . , M is known; in the step (3.2.1), using state transition priori is the function, q(x_k|x_k-1, y_0:k)=p(x_k|x_k-1), approximating the distribution of x_k by the particles, p(x_k-1|y_0:k-1)≅Σ_i=1 ^{N s} ω_k-1 ⁱδx^k _i, and obtaining:

p _m(y _k |y _0:k-1)≈Σ_i=1 ^{N s} ω_k-1 ⁱ p _m(y _k |x _k ⁱ)

where in the step (4), adopting a particle filter algorithm for the state estimation, and calculating p(x_k|H_k=m, y_0:k) and p(H_k=m|y_0:k) based on particles.

8. The computer-implemented adaptive brain-computer interface decoding method based on multi-model dynamic ensemble and executed by a processor according to claim 4, wherein α=0.9.

Images