Multi-day retraining-free motion decoding method based on muscle-bone model under variable load

Abstract

The invention discloses a multi-day retraining-free motion decoding method based on a muscle-bone model under a variable load, and relates to the technical field of biomedical engineering. The method comprises the following steps: step 1, construction of a state equation; step 2, construction of an observation equation; step 3, offline experiment; step 4, offline model training; according to the method, the state prediction error of the muscle-bone model can be corrected, the interday stability of the muscle collaborative decomposition result is improved, and the decoding performance of the muscle-bone model under the conditions of multiple days and variable loads is further improved.

Description

Technical Field

[0001] This invention relates to the field of biomedical engineering technology, and in particular to a multi-day retraining-free motion decoding method based on a musculoskeletal model under variable load. Background Technology

[0002] In the field of electromyography (EMG) interfaces and EMG prosthesis control, improving the robustness of decoding algorithms is crucial for enhancing prosthesis control performance. During daily use of prostheses, the stability of EMG signals is affected by various factors, including changes in signal characteristics over prolonged wear and variations in load weight. As a time-varying and non-stationary electrophysiological signal, these factors can alter muscle activation patterns, affecting the characteristic distribution of EMG signals. Ultimately, this leads to a decline in the decoding performance of the EMG interface, impacting the control effectiveness of the EMG prosthesis and the user experience. In daily life, prolonged prosthesis wear and variations in load weight are unavoidable; therefore, EMG interfaces need to maintain reliable decoding performance in complex and changing scenarios. This is of great significance for improving the practicality of EMG prostheses and enhancing the quality of life for amputees.

[0003] To mitigate the impact of these two factors on the performance of electromyography (EMG) interfaces, the academic community has conducted extensive research, primarily focusing on the following aspects: First, expanding the training dataset to cover diverse environmental conditions and motion tasks to increase the model's adaptability; second, introducing online learning or incremental learning strategies to enable EMG interfaces to gradually adapt to new signal characteristics during actual use; and third, utilizing dynamic neural network models, such as long short-term memory networks and recurrent neural networks, combined with time series analysis to capture the dynamic characteristics of EMG signals changing over time, providing better predictive performance and adaptability. Compared to data-driven methods that rely on large amounts of training data and complex machine learning algorithms, musculoskeletal models exhibit strong decoding robustness in application scenarios involving long-term use of EMG interfaces and varying load weights, helping to reduce the training and usage burden on users. This robustness mainly stems from the fact that musculoskeletal models are constructed based on the physiological structure and biomechanical characteristics of the human neuromuscular-skeleton system, enabling them to explain muscle contraction, stretching, and other mechanical behaviors, as well as the neuromuscular control mechanism, from a physiological perspective. Simultaneously, the biomechanical characteristics within the model are characterized by a relatively stable biomechanical model, and the model parameters have clear physical meanings, such as muscle fiber length and elastic coefficient. These parameters can only vary within a reasonable physiological range, thus ensuring the stability of the model in the face of changes in the external environment. In addition, the musculoskeletal model also considers the dynamic interaction between multi-degree-of-freedom joint movements, taking into account both the coordination of multi-joint movements and strictly following the basic principles of physics, ensuring that the entire decoding process conforms to both mechanical equilibrium and physiological requirements.

[0004] However, current research on musculoskeletal models primarily validates the inherent robustness of traditional electromyography (EMG) amplitude-based models in practical applications. While these studies demonstrate the stability of musculoskeletal models under multi-day use and variable load conditions, they do not delve into the impact of interfering factors on muscle activation patterns and musculoskeletal system biomechanical properties as deeply as data-driven approaches. Furthermore, existing research generally employs an "open-loop" prediction framework, lacking error correction mechanisms, resulting in uncorrected prediction errors and thus limiting the decoding performance of musculoskeletal models in variable application scenarios. Therefore, to further improve the decoding accuracy and robustness of musculoskeletal models while reducing the user's training burden, the following research is still needed: in-depth investigation of the dynamic changes in muscle activation and musculoskeletal system biomechanical properties during multi-day, variable load applications of prostheses, and development of new algorithms and models to adapt to these changes. Summary of the Invention

[0005] In view of the above-mentioned deficiencies of the prior art, the technical problem to be solved by the present invention is how to improve the decoding performance of the musculoskeletal model under heavy load and cross-day similar motion tasks.

[0006] To achieve the above objectives, this invention provides a multi-day retraining-free motion decoding method based on a musculoskeletal model under variable load, comprising the following steps:

[0007] Step 1: Construction of the state equation

[0008] The state equation is constructed based on a variable stiffness musculoskeletal model driven by co-activation information, with the activation level a of muscle j at time k as the input. j (k), the output state variables, i.e., the joint kinematic parameters, are...

[0009]

[0010] The first three and the last three parameters are the flexion and extension angles, angular velocity, and angular acceleration of the metacarpophalangeal joint and the wrist joint, respectively.

[0011] Step 2: Construction of observation equations

[0012] The mapping relationship between state variables and electromyographic features is established based on the back propagation neural network (BPNN), and the measurement equation is constructed. The set of electromyographic features that are robust to load changes is used as the observation value of the observation equation.

[0013] Step 3, Offline Experiment

[0014] Through the first day of the multi-day experiment, the surface electromyography signals and joint angles during the wrist and hand flexion and extension movements with two degrees of freedom were collected simultaneously, and the model training was completed.

[0015] Step 4: Offline model training

[0016] The offline model training dataset consists of three independent trials under zero-load conditions on the first day of the experiment, with each trial corresponding to a motion task in the experimental paradigm.

[0017] Step 5: Online decoding.

[0018] Furthermore, in step 1, joint stiffness that varies with muscle activation is introduced into the musculoskeletal model of the state equation, and the expression of the state equation is:

[0019]

[0020] Among them, the muscle activation level a of each muscle j (k) As input, typically, the amplitude of the electromyographic signal is used to characterize muscle activation in the musculoskeletal model. j The nonnegative matrix factorization algorithm NMF-CL2WH with added sparsity constraints was used to calculate muscle activation based on a cooperative mechanism from multi-channel surface electromyography signals. Then, the muscle activation was used to calculate the torque generated by the muscle on the d-th joint (metacarpophalangeal joint: d=1, wrist joint: d=2). and joint stiffness It is the joint damping coefficient, I d It is the moment of inertia, T s =0.01 seconds is the sampling period, w(k) is Gaussian white noise, and its covariance matrix is ​​Q.

[0021] Further, step 1 includes the following steps:

[0022] Step 1.1 Muscle activation extraction based on NMF-CL2WH

[0023] The NMF-CL2W algorithm transforms the physiological principles of muscle synergy and the muscle activation characteristics of multi-degree-of-freedom tasks into mathematical constraints, which are directly applied to the solution process of NMF. Based on the traditional NMF, the NMF-CL2W algorithm introduces sparsity constraints and L2 norm regularization terms for the co-activation weight matrix W and the co-activation coefficient matrix H, and introduces the Hadamard product to apply sparsity constraints to specified elements in matrices W and H.

[0024] Step 1.2 Construction of the variable stiffness musculoskeletal model

[0025] Further, step 1.2 includes the following steps:

[0026] Step 1.2.1 Calculate the torque of muscle pulling on the joint based on Hill's muscle mechanics model:

[0027]

[0028] Among them, the magnitude of the contractile force of muscle j is

[0029]

[0030] It is the maximum voluntary contractile force of the muscle, l j It is the current muscle fiber length, v j f represents the current contraction velocity of the muscle fiber. a (l j ), f(v j ), f p (l j Define the relationships between active contractile force and muscle fiber length, active contractile force and muscle fiber contraction velocity, and passive contractile force and muscle fiber length, respectively. j,d Let be the lever arm length of muscle j about the axis of rotation of joint d;

[0031] Step 1.2.2 Combining muscle stiffness Calculate joint stiffness:

[0032]

[0033] Muscle stiffness is equivalent to the parallel stiffness of the muscle's active elastic element (CE) and parallel elastic element (PE), and the calculation formula is as follows:

[0034]

[0035] Wherein, the stiffness of CE is

[0036]

[0037] The stiffness of PE is

[0038] It is the optimal muscle fiber length.

[0039] Further, step 2 includes the following steps:

[0040] Step 2.1 Selection of Observations

[0041] The difference between the first two principal components of the electromyography feature set that is robust to changes in load weight is used as the observed value of the observation equation. The electromyography feature set includes four electromyography time-domain features: sample entropy, zero crossings, waveform slope sign changes, and electromyography variance.

[0042] Step 2.2 Construction of observation equations

[0043] z(k) = h(x(k)) + v(k),

[0044] Where v(k) is Gaussian white noise with a covariance matrix of R, and the function h(x(k)) in the observation equation is a BPNN-based mapping model to establish the relationship between kinematic parameters and electromyographic features; a three-layer BPNN is used to construct a mapping model between joint motion parameters and electromyographic features, and this model is used as the observation equation, wherein the joint motion parameters are joint angle and angular velocity;

[0045] The training process of BPNN is as follows: using the training data collected under zero-load experiment, with joint kinematic data x as input and electromyographic features z as output, the BPNN in the observation equation is trained.

[0046] The BPNN has a three-layer structure: an input layer, a hidden layer, and an output layer. The hidden layer uses the hyperbolic tangent sigmoid function (tansig function) and has 10 neurons. The output layer uses the linear function purelin.

[0047] Further, step 2.1 includes the following steps:

[0048] Step 2.1.1 The preprocessed electromyographic signal is divided into a series of rectangular windows with a window length of 100 milliseconds and a sliding step of 10 milliseconds, and the corresponding electromyographic temporal features are calculated in each time window.

[0049] Step 2.1.2 Decentralize the electromyographic time-domain feature data. The decentralization parameters include column mean and principal component coefficients. The decentralization parameters are calculated based on the training data and are used for model validation.

[0050] Step 2.1.3 Perform principal component analysis (PCA) on the extracted electromyography time-domain feature set, and select the difference between the first two principal components as the observed value z(k).

[0051] Further, in step 3, the setup and paradigm of the offline experiment are as follows: the subject keeps their upper arm hanging naturally, elbow bent at 90 degrees, hand relaxed, and the initial state is that the palm is perpendicular to the ground; the offline experiment includes wrist flexion and extension and metacarpophalangeal joint flexion and extension; the offline experiment includes task 1, task 2 and task 3, where task 1 is metacarpophalangeal joint flexion and extension alone, task 2 is wrist flexion and extension alone, and task 3 is simultaneous flexion and extension of the metacarpophalangeal joint and wrist joint; under no external load conditions, the subject completes task 1, task 2 and task 3 at a self-selected speed and force; each task is repeated three times, and each test lasts for 20 seconds; after the first day of the experiment, an electromyographic signal under maximum voluntary contraction is collected, and after the offline experiment data is collected, the electromyographic electrode positions are marked on the forearm skin with a marker.

[0052] Furthermore, step 4 includes the following steps:

[0053] Step 4.1 Optimization of muscle physiological parameters in the musculoskeletal model

[0054] In the state equation of step 1, a constrained global optimization method is used to determine some muscle parameters, which are the optimal muscle fiber lengths. Maximum spontaneous contraction force Lever arm length r j,d Based on the model training data collected under zero-load conditions on the first day of the experiment, with the optimization objective of minimizing the sum of squared errors between the actual measured angle and the model predicted angle, the subject-specific muscle parameters and state equation expressions were determined. The model training data consisted of electromyographic signals and actual joint angles.

[0055] Step 4.2 Training the BPNN in the observation equation

[0056] Using joint kinematic data as input and electromyographic features z as output, the BPNN in the observation equation is trained; Step 4.3 Determination of unscented Kalman filter parameters.

[0057] α controls the dispersion of the Sigma points, α = 0.05; β is a non-negative constant used to describe the distribution information of the state variables, β = 0.1, and the initial state variables are set to 0. 6×1 The initial covariance matrix, the process noise covariance matrix Q, and the observation noise covariance matrix R are fully adjusted to estimate the effectiveness of different subjects' actual states.

[0058] Furthermore, in step 5, seven subjects participated in an online TAC test experiment at the same time on the first, fourth, and seventh days of the seven-day period. The online test duration was approximately 2.5 hours per day. During the daily test experiment, the subjects were required to complete the test experiment under four different external loads. A GUI was provided on the computer screen in front of the subjects, providing real-time feedback on the joint angles predicted by the decoding algorithm. The subjects controlled the virtual hand to repeatedly move from the initial position to each red target position three times by activating wrist flexion and extension and metacarpophalangeal joint flexion and extension. The four different external loads were W1 = 0 kg, W2 = 0.5 kg, W3 = 0.5 kg, W4 = 0.5 kg, W5 = 0.5 kg, W6 = 0.5 kg, W7 = 0.5 kg, W8 = 0.5 kg, W9 = 0.5 kg, W1 = 0.5 kg, W2 ...1 = 0.5 kg, W2 = 0.5 kg, W1 = 0.5 kg, W2 = 0.5 kg, W1 = 0.5 =1.0kg, W4=1.5kg. The online TAC test is divided into four parts, each corresponding to a specific load weight. The real-time control effect of each subject on the virtual hand is quantitatively evaluated through three online indicators: completion time, path efficiency, and completion rate. The completion time is the time to successfully control the virtual hand to reach the target. The path efficiency is the ratio of the shortest path from the initial position to the target position to the actual path length in the Cartesian space of joint angles. The completion rate is the ratio of the number of successful trials to the total number of trials. Both path efficiency and completion rate are dimensionless statistics.

[0059] Furthermore, the aforementioned multi-day no-retraining motion decoding method based on a musculoskeletal model under variable load can be applied to the fields of intelligent prostheses, rehabilitation equipment, and robot control.

[0060] Compared with the prior art, the present invention has the following advantages: it can correct the prediction error of the musculoskeletal model state, significantly improve the inter-day stability of the co-element decomposition results, better describe the contractile mechanical properties of the musculoskeletal system under variable load conditions, and improve the stability of decoding performance under multi-day and variable load conditions.

[0061] The following will further explain the concept, specific structure, and technical effects of the present invention in conjunction with the accompanying drawings, so as to fully understand the purpose, features, and effects of the present invention. Attached Figure Description

[0062] Figure 1 This is a schematic diagram of a multi-day retraining-free motion decoding method based on a musculoskeletal model under variable load, according to an embodiment of the present invention.

[0063] Figure 2 The offline experimental setup for this embodiment of the invention includes the positions of electromyography electrodes, angle sensors, and weighted wristbands.

[0064] Figure 3 The online TAC experiment and visual feedback for subjects in this embodiment of the invention:

[0065] (a) Online experimental scenario for subjects,

[0066] (b) Visual feedback on a screen in front of the healthy limb subject.

[0067] (c) Visual feedback on the screen in front of the healthy limb subject.

[0068] (d) Target location of the healthy limb subject during the online task;

[0069] Figure 4 The online decoding performance of three methods under different experimental days and load weight conditions in the embodiments of the present invention is shown. Detailed Implementation

[0070] The preferred embodiments of the present invention are described below with reference to the accompanying drawings to make the technical content clearer and easier to understand. The present invention can be embodied in many different forms, and the scope of protection of the present invention is not limited to the embodiments mentioned herein.

[0071] In the accompanying drawings, components with the same structure are indicated by the same numerical designation, and components with similar structures or functions are indicated by similar numerical designations. The dimensions and thicknesses of each component shown in the drawings are arbitrary, and the present invention does not limit the dimensions and thicknesses of each component. To make the illustrations clearer, the thickness of some components has been appropriately exaggerated in the drawings.

[0072] This invention addresses the problem that existing musculoskeletal models employ an "input-output" open-loop prediction framework, which cannot effectively correct motion prediction errors and suffers from severe performance degradation under long-term use or heavy load conditions. It introduces the concept of mutual correction between model predictions and measured data from Kalman filtering, constructing observation equations through a backpropagation neural network to establish a mapping relationship between state variables (joint kinematic parameters) and robust electromyographic (EMG) features to load changes, aiming to correct state prediction errors. Furthermore, existing musculoskeletal models for motion decoding often use EMG amplitude to measure muscle activation; however, this feature often lacks consistency when performing similar motion tasks across days, affecting the long-term decoding performance of musculoskeletal models. An improved nonnegative matrix factorization algorithm (NMF-CL2WH) is proposed. By effectively extracting long-term stable cofactors, it provides more stable muscle activation information for musculoskeletal models, aiming to achieve multi-day retraining-free decoding of musculoskeletal models. To address the issue that the changes in joint stiffness caused by the co-contraction of antagonistic muscle groups under high load are often overlooked, which limits the decoding performance of musculoskeletal models under dynamic load conditions, this paper proposes to introduce the variable stiffness characteristics of joints into the contractile mechanics model in order to effectively address the changes in muscle co-contraction levels caused by variations in load weight.

[0073] 1. The model part of the embodiments of the present invention includes state equations and measurement equations.

[0074] The state equation is constructed based on a variable stiffness musculoskeletal model driven by co-activation information (e.g.) Figure 1As shown in the A-state equation, the input is the activation level a of muscle j at time k. j (k), the output is the state variable (i.e., the joint kinematic parameters are)

[0075]

[0076] The first three and last three parameters are the flexion-extension angle, angular velocity, and angular acceleration of the metacarpophalangeal and wrist joints, respectively. The measurement equation is based on a backpropagation neural network (BPNN) to establish a mapping relationship between state variables and electromyographic features (e.g., ...). Figure 1 (As shown in the B-state equation).

[0077] 1.1 Construction of State Equations

[0078] By analyzing the impact of load changes on muscle activation and the biomechanical properties of the musculoskeletal system, joint stiffness that varies with muscle activation is introduced into the musculoskeletal model based on the equation of state. This addresses the changes in muscle co-contraction levels caused by load changes, thereby better describing the contractile mechanical properties of the musculoskeletal system under heavier loads and improving the decoding performance of the musculoskeletal model under load changes.

[0079] The expression for the state equation is:

[0080]

[0081] Among them, the muscle activation level a of each muscle at time k j (k) As input, typically, the amplitude of the electromyographic signal is used to characterize muscle activation in the musculoskeletal model. j This embodiment employs a nonnegative matrix factorization algorithm (NMF-CL2WH) with added sparsity constraints to calculate muscle activation based on a cooperative mechanism from multi-channel surface electromyography (EMG) signals. Subsequently, the muscle activation is used to calculate the torque generated by the muscle on the d-th joint (metacarpophalangeal joints: d=1, wrist joints: d=2). and joint stiffness It is the joint damping coefficient, I d It is the moment of inertia, T s =0.01 seconds is the sampling period, w(k) is Gaussian white noise, and its covariance matrix is ​​Q.

[0082] 1) Muscle activation extraction based on NMF-CL2WH

[0083] This invention proposes to use the NMF-CL2W algorithm to transform the physiological principles of muscle synergy and the muscle activation characteristics of multi-degree-of-freedom tasks into mathematical constraints, which are then directly applied to the solution process of NMF.

[0084] Building upon traditional NMF, this embodiment proposes the NMF-CL2WH algorithm. By introducing sparsity constraints and L2 norm regularization terms for the co-activation weight matrix W and the co-activation coefficient matrix H, these constraints guide the solution process to avoid solutions with small reconstruction errors but inconsistent with physiological / task logic. Instead, the algorithm favors solutions that meet both accuracy requirements and the expectation of a stable co-activation structure, thus optimizing the NMF decomposition results and improving the determinism of the decomposition. To impose motion activation-related constraints on W and H, this embodiment introduces the concept of the Hadamard product, used to apply sparsity constraints to specified elements in the matrices. For two matrices of the same order, X and Y, their Hadamard product is denoted as... If the dimensions of X and Y are both m×n, then their Hadamard product is... It is also an m×n matrix, where the element in the i-th row and j-th column is the product of the corresponding elements of X and Y. The sparsity constraint matrix is ​​defined as C. w and C h Their sizes are the same as W and H, respectively, and they are used to constrain W and H.

[0085] First, this embodiment involves four coercive elements, all satisfying a Variation Account For parameter greater than 0.9. The joint labels d for the metacarpophalangeal and wrist joints are defined as 1 and 2, respectively, and the coercive element labels for the flexion and extension movements of each joint are + and -, respectively. Therefore, the electromyographic feature data Z can be decomposed into W×H, with each column and row in W and H defined as follows:

[0086]

[0087] Where t represents time and T is the matrix transpose symbol.

[0088] a)C w Determination method

[0089] In this embodiment, a total of 6 channels of electromyography signals were collected in the experiment, and the weight of each coercive element is:

[0090]

[0091] Where j represents the j-th muscle. To compare the activation contributions of different muscles within the same coercive element, the activation weight matrix corresponding to each coercive element is normalized using the maximum weight of that coercive element:

[0092]

[0093] To compare the activation percentage of the same muscle across different coercors, P is introduced. d The activation percentage is represented by the following calculation method:

[0094] (1) Activation coefficient and Normalize each value separately; the normalized activation coefficients are as follows: and

[0095] (2) and Transformed into and Ensure and The product of these terms is still equal to Z, where... So and The calculation methods for each element are as follows:

[0096]

[0097] (3) The formula for calculating the activation percentage of muscle among synergists is as follows:

[0098]

[0099] In this embodiment, regarding the above weighted correlation coefficients A thresholding method was used to determine the activation contribution of muscles to different movements. This was determined by the following method: if and only if and If all values ​​are greater than 0.2, then the j-th muscle is considered to contribute to the d-th degree of freedom of flexion. otherwise Similarly, and If all values ​​are greater than 0.2, then the j-th muscle is considered to contribute to the extension of the d-th degree of freedom. otherwise

[0100] b)C h Determination method

[0101] constraint matrix C h The activation coefficient matrix H of the co-elements is applied. When the d-th degree of freedom participates in the motion, the element corresponding to the d-th degree of freedom of the co-element in the constraint matrix is ​​1, while the element corresponding to the co-element that does not participate in the motion is 0. Thus, based on the task type of the selected training data, C can be determined in advance. h .

[0102] c) L2 norm regularization term

[0103] To avoid the algorithm relying excessively on the electromyographic data used during calibration, and considering that the L2 norm is an important indicator affecting the algorithm's structural complexity, in order to simplify the model and improve the algorithm's robustness and the reliability of the electromyographic control method, in C... w and C h Based on the constraints, add L2 norm regularization terms for the row vectors W and H.

[0104] d) The objective function of NMF-CL2WH is

[0105]

[0106] In this context, the subscript Fro denotes the Frobenius norm, and λ0 and β0 are regularization coefficients used to weigh the impact of the calibration feature dataset on the model. These coefficients were selected through offline cross-validation analysis, with values ​​of 5 and 1. express The j-th row of the matrix, Representation matrix The lth column.

[0107] e) Muscle activation calculation

[0108] With the training dataset, subject-specific co-activation weight matrices can be obtained. And used to solve the activation coefficient matrix H of the cooperational elements. M :

[0109]

[0110] Where H M The four row vectors represent the electromyographic components that activate the metacarpophalangeal joint flexion and extension, and the wrist joint flexion and extension, respectively, and are further converted into muscle activation to drive the musculoskeletal model. In mathematics, st means "constrained by", followed by a constraint condition.

[0111] 2) Construction of variable stiffness musculoskeletal model

[0112] Calculate the torque by which muscles pull on joints based on Hill's muscle mechanics model:

[0113]

[0114] Among them, the magnitude of the contractile force of muscle j is

[0115]

[0116] It is the maximum voluntary contractile force of the muscle, l j It is the current muscle fiber length, v j This indicates the current contraction velocity of the muscle fiber.a (l j ), f(v j ), f p (l j The relationships between active contractile force and muscle fiber length, active contractile force and muscle fiber contraction velocity, and passive contractile force and muscle fiber length were defined respectively. j,d Let be the lever arm length of muscle j about the rotation axis of joint d. This should be considered in conjunction with muscle stiffness. Calculate joint stiffness:

[0117]

[0118] Muscle stiffness is equivalent to the parallel stiffness of the muscle's active elastic element (CE) and parallel elastic element (PE), and the calculation formula is as follows:

[0119]

[0120] Wherein, the stiffness of CE is

[0121]

[0122] The stiffness of PE is

[0123] It is the optimal muscle fiber length.

[0124] 1.2 Construction of Observation Equations

[0125] To minimize errors caused by modeling uncertainties and improve the motion estimation accuracy of the musculoskeletal model under multi-day and external load weight variations, this embodiment constructs an observation equation for the variable stiffness musculoskeletal model. The electromyographic feature set robust to load changes is used as the observation value of the observation equation. The observation equation establishes a mapping relationship between state variables and electromyographic features using an artificial neural network. Through an unscented Kalman filter algorithm, the state prediction value of the musculoskeletal model is optimally fused with the sensor observation data to correct the prediction error of the "open-loop" state equation.

[0126] 1) Selection of observation values

[0127] During motion decoding, the input signal consists only of electromyographic (EMG) signals acquired by surface EMG electrodes; no other sensors can provide information on joint movement and load weight. This embodiment proposes an EMG feature set robust to load weight changes, using the difference between its first two principal components as the observation value. This feature set consists of four EMG time-domain features: sample entropy, zero crossings, waveform slope sign changes, and EMG variance. To extract these features, the preprocessed EMG signal is segmented into a series of rectangular windows with a window length of 100 ms and a sliding step of 10 ms, and the corresponding EMG features are calculated within each time window. Subsequently, to reduce the complexity of high-dimensional features and retain important information, principal component analysis (PCA) is performed on the extracted time-domain feature set, and the difference between the first two principal components is selected as the observation value z(k). It is important to note that when performing PCA analysis, the data first needs to be decentralized. The decentralization parameters (column mean and principal component coefficients) are calculated based on the training data and used for model validation.

[0128] 2) Construction of observation equations

[0129] z(k) = h(x(k)) + v(k),

[0130] Where v(k) is Gaussian white noise with a covariance matrix of R. The function h(x(k)) in the observation equation is a BPNN-based mapping model that establishes the relationship between kinematic parameters and electromyographic features. A three-layer BPNN is used to construct a mapping model between joint motion parameters (joint angles and angular velocities) and electromyographic features, and this model is used as the observation equation. The training process of the BPNN is as follows: using the training data collected under zero-load experiment on the first day of the experiment, with joint kinematic data x as input and electromyographic features z as output, the BPNN in the observation equation is trained. The BPNN has a three-layer structure: an input layer, a hidden layer, and an output layer. The hidden layer uses the hyperbolic tangent sigmoid function (tansig function), and the number of neurons in the hidden layer is 10. The output layer uses the linear function purelin.

[0131] 2. The experimental part of this invention includes offline experiments, offline model training, and online decoding.

[0132] 2.1 Offline Experiment

[0133] like Figure 2 As shown, through offline experiments under zero-load conditions on the first day of a multi-day experiment, surface electromyographic signals and joint angles during wrist and hand flexion and extension movements of two degrees of freedom were collected simultaneously, and model training was completed.

[0134] The specific experimental setup and paradigm were as follows: Subjects kept their upper arms hanging naturally with elbows bent at 90 degrees and hands relaxed, initially with their palms perpendicular to the ground. The experimental movements included wrist flexion and extension, and metacarpophalangeal joint flexion and extension. Specific tasks were as follows: Task 1 was metacarpophalangeal joint flexion and extension alone; Task 2 was wrist flexion and extension alone; and Task 3 was simultaneous flexion and extension of both the metacarpophalangeal and wrist joints. Under no external load, subjects were required to complete the three experimental movements at their chosen speed and force. Each movement was repeated three times, with each trial lasting 20 seconds. At the end of the first day of the experiment, electromyography (EMG) signals were collected under maximal voluntary contraction for EMG signal normalization. After offline training data collection, the EMG electrode positions were marked on the forearm skin using a marker to minimize the impact of electrode misalignment on subsequent experiments.

[0135] 2.2 Offline Model Training:

[0136] Based on the aforementioned offline experimental setup and paradigm, this embodiment proposes a decoding method based on a musculoskeletal model. The training dataset for this model consists of three independent trials under zero-load conditions on the first day of the experiment, with each trial corresponding to a different motion task within the experimental paradigm. The specific model training mainly includes three parts:

[0137] 1) Optimization of muscle physiological parameters in musculoskeletal model

[0138] In the above state equation, since some muscle parameters cannot be directly measured, a constrained global optimization method was used to determine these parameters, including the optimal muscle fiber length. Maximum spontaneous contraction force Lever arm length r j,d By combining model training data (electromyographic signals and actual joint angles) collected under zero-load experimental conditions, and with the optimization objective of minimizing the sum of squared errors between the actual measured angles and the model-predicted angles, the subject-specific muscle parameters and state equation expressions were finally determined.

[0139] 2) Training of BPNN in the observation equation

[0140] Using joint kinematics data as input and electromyographic features z as output, the BPNN in the observation equation is trained.

[0141] 3) Determination of Unscented Kalman Filter Parameters

[0142] α controls the dispersion of the Sigma points, and is set to 0.05. β is a non-negative constant used to describe the distribution information of the state variables, and is set to 0.1. The initial state variables are set to 0. 6×1The initial covariance matrix, the process noise covariance matrix Q, and the observation noise covariance matrix R need to be fully adjusted based on the estimation results of different subjects' actual states.

[0143] 2.3 Online Decoding

[0144] After completing offline model training, seven subjects participated in online TAC tests at the same time on days one, four, and seven within a seven-day period (approximately 2.5 hours of online experiment time per day) to simulate a multi-day real-time use scenario for the electromyography interface. During each day's experiment, subjects were required to complete the experiment under different external loads to comprehensively evaluate the decoding performance of different decoding methods in a multi-day, no-retraining scenario under varying loads. During this process, a GUI was provided on the computer screen in front of the subjects, providing real-time feedback on the joint angles predicted by the decoding algorithm, such as... Figure 3 As shown. Subjects need to control the virtual hand to repeatedly move from the initial position to each red target position three times by activating the wrist flexion and extension and metacarpophalangeal joint flexion and extension. Healthy limb subjects need to complete the online TAC experiment under four different external load conditions every day: W1=0kg, W2=0.5kg, W3=1.0kg, W4=1.5kg. The daily online TAC experiment is divided into four parts, each corresponding to a specific load weight. Under each load weight, the experiment is further divided into three groups, each corresponding to a decoding method: namely the method proposed in this embodiment, linear regression, and the musculoskeletal model based on electromyography amplitude (CROUCH DL, HUANG H. Lumped-parameter electromyogram-driven musculoskeletal hand model: a potential platform for real-time prosthesis control[J]. Journal of Biomechanics,2016,49(16):3901-3907.). Therefore, subjects need to complete 4*3=12 sets of experiments every day, and each set of experiments contains five targets. The real-time control performance of each subject on the virtual hand is quantitatively evaluated using three online metrics: completion time (the time to successfully control the virtual hand to reach the target), path efficiency (the ratio of the shortest path length from the initial position to the target position to the actual path length in Cartesian space of joint angles), and completion rate (the ratio of the number of successful trials to the total number of trials). These three evaluation metrics are widely used to assess the online control performance of electromyographic decoding algorithms. Path efficiency and completion rate are dimensionless statistics; the closer the value is to 100% and the shorter the completion time, the stronger the subject's real-time control ability on the virtual hand, i.e., the better the control performance of the electromyographic interface.

[0145] 3. Evaluation of the online decoding performance of the embodiments of the present invention

[0146] All model training for the comparison methods was completed on the first day of the experiment. Subsequent online control experiments directly used the models trained on the first day, without retraining for different experiment dates or load weights. Figure 4 The table lists the online decoding performance of three methods under different experimental days and load conditions, including this embodiment, the musculoskeletal model based on electromyography amplitude (General-EMGMM), and linear regression (LR). For commonly used degrees of freedom in electromyographic prosthetic control (hand grasping and wrist flexion / extension), the online TAC experiments conducted in this embodiment under multiple days and different load weights show that the proposed method can complete offline model training using the first day's no-load experimental data. It can achieve stable continuous wrist and hand motion decoding without retraining under subsequent experimental conditions, and its decoding performance is not easily affected by load changes or cross-day applications. Compared with the musculoskeletal model based on electromyography amplitude and the linear regression method, the proposed method performs best in key indicators such as completion time, completion rate, and path efficiency, demonstrating superior decoding accuracy and robustness. This fully proves the effectiveness and superiority of the proposed method in multi-day no-retraining scenarios under varying loads.

[0147] The method proposed in this invention only requires initial model training by collecting training data under no additional load during the first use. No additional training is required in subsequent multi-day and variable load scenarios. It can achieve robust decoding of continuous flexion and extension movements of the wrist and hand with two degrees of freedom across days (one week) and across loads (0-1.5kg load), thereby significantly reducing the training burden on subjects. These breakthroughs will promote the practical application of myoelectric interfaces based on musculoskeletal models and provide a more stable and reliable control solution for myoelectric prostheses and other devices that require myoelectric auxiliary control.

[0148] The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and variations based on the concept of the present invention without creative effort. Therefore, all technical solutions that can be obtained by those skilled in the art based on the concept of the present invention through logical analysis, reasoning, or limited experimentation on the basis of existing technology should be within the scope of protection defined by the claims.

Claims

1. A multi-day no-retraining motion decoding method based on a musculoskeletal model under variable load, characterized in that, Includes the following steps: Step 1: Construction of the state equation The state equation is constructed based on a variable stiffness musculoskeletal model driven by co-activation information, with the activation level a of muscle j at time k as the input. j (k), the output state variables, i.e., the joint kinematic parameters, are... The first three and the last three parameters are the flexion and extension angles, angular velocity, and angular acceleration of the metacarpophalangeal joint and the wrist joint, respectively. Step 2: Construction of observation equations The mapping relationship between state variables and electromyographic features is established based on the back propagation neural network (BPNN), and the measurement equation is constructed. The set of electromyographic features that are robust to load changes is used as the observation value of the observation equation. Step 3, Offline Experiment Through the first day of the multi-day experiment, the offline experiment under zero-load conditions was conducted to simultaneously collect surface electromyographic signals and joint angles during the two-degree-of-freedom flexion and extension movements of the wrist and hand, and to complete the model training. Step 4: Offline model training The offline model training dataset consists of three independent trials under zero-load conditions on the first day of the experiment, with each trial corresponding to a motion task in the experimental paradigm. Step 5: Online decoding.

2. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 1, characterized in that, In step 1, joint stiffness that varies with muscle activation is introduced into the musculoskeletal model of the state equation. The expression for the state equation is: Among them, the muscle activation level a of each muscle j (k) As input, typically, the amplitude of the electromyographic signal is used to characterize muscle activation in the musculoskeletal model. j The nonnegative matrix factorization algorithm NMF-CL2WH with added sparsity constraints was used to calculate muscle activation based on a cooperative mechanism from multi-channel surface electromyography signals. Then, the muscle activation was used to calculate the torque generated by the muscle on the d-th joint (metacarpophalangeal joint: d=1, wrist joint: d=2). and joint stiffness It is the joint damping coefficient, I d It is the moment of inertia, T s =0.01 seconds is the sampling period, w(k) is Gaussian white noise, and its covariance matrix is ​​Q.

3. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 2, characterized in that... Step 1 includes the following steps: Step 1.1 Muscle activation extraction based on NMF-CL2WH The NMF-CL2W algorithm transforms the physiological principles of muscle synergy and the muscle activation characteristics of multi-degree-of-freedom tasks into mathematical constraints, which are directly applied to the solution process of NMF. Based on the traditional NMF, the NMF-CL2W algorithm introduces sparsity constraints and L2 norm regularization terms for the co-activation weight matrix W and the co-activation coefficient matrix H, and introduces the Hadamard product to apply sparsity constraints to specified elements in matrices W and H. Step 1.2 Construction of the variable stiffness musculoskeletal model.

4. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 3, characterized in that, Step 1.2 includes the following steps: Step 1.2.1 Calculate the torque of muscle pulling on the joint based on Hill's muscle mechanics model: Among them, the magnitude of the contractile force of muscle j is It is the maximum voluntary contractile force of the muscle, l j It is the current muscle fiber length, v j f represents the current contraction velocity of the muscle fiber. a (l j ), f(v j ), f p (l j Define the relationships between active contractile force and muscle fiber length, active contractile force and muscle fiber contraction velocity, and passive contractile force and muscle fiber length, respectively. j,d Let be the lever arm length of muscle j about the axis of rotation of joint d; Step 1.2.2 Combining muscle stiffness Calculate joint stiffness: Muscle stiffness is equivalent to the parallel stiffness of the muscle's active elastic element (CE) and parallel elastic element (PE), and the calculation formula is as follows: Wherein, the stiffness of CE is The stiffness of PE is It is the optimal muscle fiber length.

5. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 1, characterized in that... Step 2 includes the following steps: Step 2.1 Selection of Observations The difference between the first two principal components of the electromyography feature set that is robust to changes in load weight is used as the observed value of the observation equation. The electromyography feature set includes four electromyography time-domain features: sample entropy, zero crossings, waveform slope sign changes, and electromyography variance. Step 2.2 Construction of observation equations; z(k) = h(x(k)) + v(k), Where v(k) is Gaussian white noise with a covariance matrix of R, and the function h(x(k)) in the observation equation is a BPNN-based mapping model to establish the relationship between kinematic parameters and electromyographic features; a three-layer BPNN is used to construct a mapping model between joint motion parameters and electromyographic features, and this model is used as the observation equation, wherein the joint motion parameters are joint angle and angular velocity; The training process of BPNN is as follows: combining the training data collected under the zero-load experiment on the first day of the experiment, with joint kinematic data x as input and electromyographic features z as output, the BPNN in the observation equation is trained. The BPNN has a three-layer structure: an input layer, a hidden layer, and an output layer. The hidden layer uses the hyperbolic tangent sigmoid function (tansig function) and has 10 neurons. The output layer uses the linear function purelin.

6. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 5, characterized in that, Step 2.1 includes the following steps: Step 2.1.1 The preprocessed electromyographic signal is divided into a series of rectangular windows with a window length of 100 milliseconds and a sliding step of 10 milliseconds, and the corresponding electromyographic temporal features are calculated in each time window. Step 2.1.2 Decentralize the electromyographic time-domain feature data. The decentralization parameters include column mean and principal component coefficients. The decentralization parameters are calculated based on the training data and are used for model validation. Step 2.1.3 Perform principal component analysis (PCA) on the extracted electromyography time-domain feature set, and select the difference between the first two principal components as the observed value z(k).

7. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 1, characterized in that, In step 3, the setup and paradigm of the offline experiment are as follows: the subject keeps their upper arm hanging naturally with the elbow bent at 90 degrees and the hand relaxed, initially with the palm perpendicular to the ground; the offline experiment includes wrist flexion and extension and metacarpophalangeal joint flexion and extension; the offline experiment includes Task 1, Task 2 and Task 3, where Task 1 is metacarpophalangeal joint flexion and extension alone, Task 2 is wrist flexion and extension alone, and Task 3 is simultaneous flexion and extension of the metacarpophalangeal joint and wrist joint; under no external load, the subject completes Task 1, Task 2 and Task 3 at a self-selected speed and force; each task is repeated three times, and each test lasts 20 seconds; after the first day of the experiment, an electromyographic signal under maximum voluntary contraction is collected, and after the offline experiment data is collected, the electromyographic electrode positions are marked on the forearm skin with a marker.

8. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 1, characterized in that, Step 4 includes the following steps: Step 4.1 Optimization of muscle physiological parameters in the musculoskeletal model In the state equation of step 1, a constrained global optimization method is used to determine some muscle parameters, which are the optimal muscle fiber lengths. Maximum spontaneous contraction force Lever arm length r j,d Based on the model training data collected under zero-load experimental conditions, with the optimization objective of minimizing the sum of squared errors between the actual measured angle and the model predicted angle, the subject-specific muscle parameters and state equation expressions are determined. The model training data consists of electromyographic signals and actual joint angles. Step 4.2 Training the BPNN in the observation equation Using joint kinematics data as input and electromyographic features z as output, the BPNN in the observation equation is trained. Step 4.3 Determining the parameters of the unscented Kalman filter α controls the dispersion of the Sigma points, α = 0.05; β is a non-negative constant used to describe the distribution information of the state variables, β = 0.1, and the initial state variables are set to 0. 6×1 The initial covariance matrix, the process noise covariance matrix Q, and the observation noise covariance matrix R are fully adjusted to estimate the effectiveness of different subjects' actual states.

9. The method for multi-day retraining-free motion decoding based on a musculoskeletal model under variable load as described in claim 1, characterized in that, In step 5, seven subjects participated in an online TAC test experiment at the same time on the first, fourth, and seventh days of the seven-day period. The online test duration was approximately 2.5 hours per day. During the daily test experiment, subjects were required to complete the test under four different external loads. A GUI was provided on the computer screen in front of the subjects, providing real-time feedback on the joint angles predicted by the decoding algorithm. Subjects controlled the virtual hand to repeatedly move from the initial position to each red target position three times by activating wrist flexion and extension and metacarpophalangeal joint flexion and extension. The four different external loads were W1 = 0 kg, W2 = 0.5 kg, W3 = 1 kg, W4 = 0 kg, W5 = 0 kg, W6 = 0 kg, W7 = 0 kg, W8 = 0 kg, W9 = 0 kg, W1 = 0 kg, W2 = 0 kg, W3 = 0 ... The online TAC test, consisting of 0kg and W4=1.5kg, is divided into four parts, each corresponding to a specific load weight. The real-time control effect of each subject on the virtual hand is quantitatively evaluated through three online indicators: completion time, path efficiency, and completion rate. The completion time is the time to successfully control the virtual hand to reach the target. The path efficiency is the ratio of the shortest path from the initial position to the target position to the actual path length in the Cartesian space of joint angles. The completion rate is the ratio of the number of successful trials to the total number of trials. Both path efficiency and completion rate are dimensionless statistics.

10. The multi-day no-retraining motion decoding method based on a musculoskeletal model under variable load as described in any one of claims 1-9, is applied in the fields of intelligent prostheses, rehabilitation equipment, and robot control.

Images