Joint Angle Prediction Method Based on Motor Units and Neuromusculoskeletal Models

Abstract

The present invention discloses a joint angle prediction method based on motor units and neuromusculoskeletal models, including: 1. Acquisition and preprocessing of electromyography data and motion data during joint movement; 2. Extracting the activation positions of motor units and clustering them using the non-negative matrix factorization algorithm; 3. Estimating neural activation through a regression model and a twitch force model; 4. Completing the calculation from neural activation to joint angle based on a specific musculoskeletal model; 5. Optimizing the extended parameter vector in the above-mentioned musculoskeletal model based on a global heuristic search algorithm, and using the optimized neuromusculoskeletal model for joint angle prediction. The present invention can decode motion behaviors from underlying neural commands according to the laws of exercise physiology, thereby enabling more robust and accurate joint angle prediction.

Description

Technical Field

[0001] The present invention belongs to the technical field of electromyogram signal processing, and particularly relates to a joint angle prediction method based on motor units and neuromusculoskeletal models, which is mainly applied to intelligent and physiological neural command and movement behavior decoding. Background Art

[0002] Limb movement is ubiquitous in life, and it is generated by the contraction of corresponding muscle fibers under the control of neural commands, which drives the movement of bones. When muscle contraction causes joint movement behavior, electromyogram signals are also generated. In essence, they are the result of the superposition of the activities of motor units (MUs) recruited by the central nervous system on the skin surface, reflecting the intention and behavior of movement. In recent years, surface electromyography (SEMG) has been widely used due to its portability and non-invasiveness. In particular, using SEMG for joint angle prediction has good application prospects in fields such as rehabilitation robots and motion control, and has received extensive attention. Currently, the joint angle prediction methods based on SEMG are mainly divided into two methods: model-free and model-based. The model-free methods mainly use machine learning, neural networks and other methods to establish the mapping relationship between SEMG data and macroscopic limb movements. However, this method relies on rich training data, and the trained model is a 'black box' with poor physiological interpretability and robustness. The model-based methods mainly perform joint angle prediction based on neuromusculoskeletal models. This method includes the physiological and anatomical structure factors of the musculoskeletal system, and can provide a clear representation between SEMG and joint movement characteristics. It has excellent robustness and physiological interpretability, and has received extensive attention in recent years. A key step in driving the neuromusculoskeletal model is to extract neural drive information from SEMG, so as to calculate subsequent muscle activities, muscle forces, joint accelerations and other information using neural drive. The current research mainly extracts the envelope of macroscopic SEMG, or uses the muscle synergy information mined by NMF as neural activation. However, this method can only obtain a rough representation of neural drive information and cannot represent the underlying neural commands more detailedly and physiologically.

[0003] With the development of high-density SEMG technology, blind source separation can be effectively used to mine the spatial information of high-density SEMG, making it possible to obtain MU activities from non-invasively measured SEMG. Currently, model-free methods based on MU activities have been widely studied and achieved excellent results. There have also been some attempts at model-based methods, such as simply superimposing all the decomposed MUs to obtain a composite firing sequence, and then using the firing rate curve obtained after low-pass filtering as the neural activation information. However, these attempts have two key problems: 1. The mapping from MUs to muscle-tendon units is not practical and not physiological. Surface electromyogram decomposition requires the rich spatial information provided by high-density SEMG. Current methods mostly use special high-density electrodes placed on specific muscle-tendon units, but this method is not practical enough. And for some highly mixed muscle groups such as the extensor and flexor muscle groups of the forearm, since the muscle-tendon units cover each other and are close in distance, it is impossible to collect the activities of only specific muscle-tendon units. Therefore, for the forearm muscle group, the collected SEMG activities are the result of the superposition of the activities of multiple muscle-tendon units, and the decomposed MU activities also come from multiple muscle-tendon units. So how to allocate the decomposed MU activities to the corresponding muscle-tendon units in a way that follows physiological laws is an urgent problem to be solved. 2. The decoding of MU activities is not sufficient by only using the composite firing information of MUs. Current model-based methods only use the firing information of MUs and do not make full use of other information. In recent years, model-free methods based on MU activities have shown that making full use of the waveform information of MUs and considering the differences between MUs has obvious advantages for the force estimation task. Therefore, how to perform more scientific and physiological decoding of MU activities in model-based methods is also an urgent problem to be solved.

[0004] For joint angle prediction based on neuromusculoskeletal models, the above two problems are the pain points that have not been solved yet. Summary of the Invention

[0005] In order to overcome the deficiencies of the existing joint angle prediction based on neuromusculoskeletal models, the present invention proposes a joint angle prediction method based on motor units and neuromusculoskeletal models, aiming to be able to decode the underlying neural commands more scientifically and physiologically, and through the musculoskeletal model, clearly represent the underlying neural commands to macroscopic movements in a way that follows physiological laws, so as to achieve more robust and accurate joint angle prediction.

[0006] In order to solve the technical problems, the present invention adopts the following technical solutions:

[0007] The joint angle prediction method based on motor units and neuromusculoskeletal models of the present invention is characterized in that it is carried out according to the following steps:

[0008] Step 1: Acquisition and preprocessing of EMG data and motion data during joint movement:

[0009] Step 1.1: Use an EMG measurement device and high-density electrodes with M channels to collect high-density surface EMG signals at time t during joint movement at the muscle to be measured, denoted as x(t) = [x1(t), x2(t), …, x m (t), …, x M (t)] T , where x m (t) represents the EMG data of the m-th channel at time t, thus obtaining a high-density surface EMG signal X with a duration of T1; T represents the transpose;

[0010] Meanwhile, use a motion acquisition device to collect the angular information of S degrees of freedom in the joint, denoted as where represents the angle of the s-th degree of freedom in the joint at time t, thus obtaining the true joint angle with a duration of T1

[0011] Use an EMG measurement device to collect the EMG peak value when the subject performs a maximum voluntary contraction, denoted as MVC;

[0012] Step 1.2: Use a blind source separation algorithm to decompose the high-density surface EMG signal X with a duration of T1, obtaining N motor unit firing sequences S = [s1, s2, …, s n , …, s N T and its two-dimensional waveform W = [w 11 , w 12 , …, w 1M , …, w nm , …, w NM T and the residual EMG waveform matrix Residual, where s n represents the n-th motor unit firing sequence with a duration of T1, and w nm represents the waveform of the n-th motor unit firing sequence s n in the m-th channel;

[0013] Step 2: Use the non-negative matrix factorization algorithm to extract the activation positions of motor units and cluster them:

[0014] Extract the activation positions of N motor units from the two-dimensional waveform W of N M channels, denoted as P = [p1, …, p n …, p N T , where p n represents the activation position of the n-th motor unit firing sequence s n ;​​​

[0015] Let the number of cluster centers be the same as the assumed number of muscle-tendon units, denoted as k. Then, use the clustering algorithm to cluster the activation positions P of N motor units to obtain k cluster centers denoted as C = [c1, …, c i , … c k T and the class indices g of each motor unit, g = [g1, …, g n …, g N T , where c i represents the two-dimensional position of the i-th cluster center, and g n represents the class index of the n-th motor unit;

[0016] Step 3: Estimate neural activation through the regression model and twitch force model:

[0017] Estimate the neural activation of k muscle-tendon units from S and W, denoted as U = [u1, u2, …, u i , …, u k T , where u i represents the neural activation of the i-th muscle-tendon unit with a duration of T1, and u i = [u i (1), u i (2), …, u i (t), …, u i (T1)], where u i (t) represents the neural activation of the i-th muscle-tendon unit at time t;

[0018] Step 4: Complete the calculation from neural activation to joint angle based on the musculoskeletal model:

[0019] Step 4.1: Nonlinearize U using Equation (1) to obtain the corresponding muscle activation z:

[0020]

[0021] In Equation (1), z is a nonlinear factor, A = [a1, a2, …, a i , …, a k T , where a i = [a i (1), a i (2), …, a i (t), …, a i (T1)], and a i (t) represents the muscle activation of the i-th muscle-tendon unit at time t; ​​​​

[0022] Step 4.2: Initialize \(t = 0\);

[0023] Set the estimated joint angle \(\theta(t)=[\theta_1(t),\theta_2(t),\ldots,\theta s (t),\ldots,\theta S (t)]\), where \(\theta s (t)\) represents the angle of the \(s\)-th degree of freedom in the estimated joint at time \(t\);

[0024] Set the muscle fiber contraction velocity \(v(t)=[v_1(t),v_2(t),\ldots,v i (t),\ldots,v k (t)]\) T as the zero vector, where \(v i (t)\) represents the muscle fiber contraction velocity of the \(i\)-th muscle-tendon unit at time \(t\);

[0025] Given the parameter vector \(h\) of the musculoskeletal model:

[0026]

[0027] In Equation (2), \(\varphi o,i , and are the maximum isometric force, optimal muscle fiber length, optimal pennation angle, tendon length, and length scaling factor of the \(i\)-th muscle-tendon unit, respectively;

[0028] Step 4.3: Calculate the muscle forces of the \(k\) muscle-tendon units at time \(t\) using the Hill muscle model where represents the muscle force of the \(i\)-th muscle-tendon unit at time \(t\);

[0029] Step 4.4: Obtain the torque \(\tau(t)\) of the \(S\) degrees of freedom in the joint at time \(t\) using Equation (3):

[0030]

[0031] In Equation (3), \(\tau(t)=[\tau_1(t),\tau_2(t),\ldots,\tau s (t),\ldots,\tau S (t)]\), \(\tau s (t)\) represents the torque of the \(s\)-th degree of freedom in the joint at time \(t\); \(r i (t)\) represents the moment arm vector of the \(i\)-th muscle-tendon unit at time \(t\), and \(r i (t)=[r i1 (t),r i2 (t),\ldots,r is (t),\ldots,r iS (t)],ris (t) is the moment arm of the i-th muscle-tendon unit in the s-th degree of freedom of the joint;

[0032] Step 4.5: Solve the state equation shown in Equation (4) using numerical integration to obtain the estimated joint angle θ(t + 1) at time t + 1;

[0033]

[0034] In Equation (4), M is the inertia matrix, C is the damping matrix, K is the stiffness matrix, and G is the gravity matrix; represents the first-order differential of the estimated joint angle at time t, represents the second-order differential of the estimated joint angle at time t;

[0035] Step 4.6: After assigning t + 1 to t, if t ≠ T1, return to Step 4.3 and execute sequentially to solve the estimated joint angle at the next moment;

[0036] Step 5: Optimize the extended parameter vector in the musculoskeletal model based on the global heuristic search algorithm:

[0037] Step 5.1: Initialize the extended parameter vector

[0038]

[0039] In Equation (5), a and b represent two regression coefficients, and c is the amplitude scaling factor;

[0040] Step 5.2: Use the extended parameter vector Process the high-density surface electromyogram signal X with a duration of T1 according to the process of Steps 2 - 4 to obtain the estimated joint angles with a duration of T1, denoted as Θ = [θ(0), θ(1), …, θ(t), …, θ(T1)];

[0041] Step 5.3: Using the true joint angle as a reference, within the set constraint range, use the global optimization algorithm to optimize and solve the extended parameter vector so that the root mean square error between the estimated joint angle Θ and the true joint angle is minimized, thereby obtaining the optimized extended parameter vector Use The determined neuromusculoskeletal model realizes the prediction of joint angles.

[0042] The characteristics of the joint angle prediction method based on motor units and neuromusculoskeletal models described in the present invention also lie in that the extraction and clustering of the motor unit activation positions in Step 2 are carried out according to the following steps:

[0043] Step 2.1. Rearrange the N two-dimensional waveforms W with M channels into N waveform matrices V1, V2, … V n , … V N , where V n represents the waveform matrix of the nth motor unit, each row of which represents a channel and each column represents a time sampling point;

[0044] Step 2.2. Decompose the waveform matrix V n of the nth motor unit by using the non-negative matrix factorization algorithm shown in Equation (6):

[0045]

[0046] In Equation (6), k is the number of muscle tendon units, Q n is the activation weight matrix of the nth motor unit, and H n is the time-varying curve matrix of the nth motor unit;

[0047] Step 2.3. For each row in H n , solve its peak value to obtain k time-varying curve peak values where δ i represents the peak value of the ith time-varying curve;

[0048] For each column in Q n , solve its peak value to obtain k activation weight peak values

[0049] Multiply the k time-varying curve peak values and the k activation weight peak values one by one in order to obtain k waveform peak values, and record the position index of the maximum value among the k waveform peak values as γ n ;

[0050] Step 2.4. Take the γ n th column of the activation weight matrix Q n as the main activation pattern, and rearrange its M activation weights into the arrangement form of the high-density electrodes;

[0051] Step 2.5. Set an activation weight threshold ε, obtain the two-dimensional positions of the M activation weights that are greater than ε and take the average to obtain the activation position p n of the nth motor unit;

[0052] Step 2.6. Obtain the activation positions p1, …, p n …, p N of all N motor units according to the process from Step 2.2 to Step 2.5;

[0053] Step 2.7: Use the clustering algorithm to cluster the activation positions P of N motor units to obtain k cluster centers denoted as C and the class indices g of each motor unit.

[0054] The said Step 3 includes:

[0055] Step 3.1: According to the cluster center positions C and the positions of the electrodes relative to the muscle tendon unit, pair the k classes with the k muscle tendon units one by one; according to the class index g of the motor units, assign the N motor unit firing sequences S to the k muscle tendon units one by one;

[0056] Step 3.2: Calculate the maximum peak values of the waveforms of N motor units in M channels from W respectively, so as to obtain N maximum peak values, and normalize the N maximum peak values using MVC to obtain N normalized maximum peak values, denoted as ρ = [ρ1, ρ2, …, ρ n , …, ρ N T , where ρ n represents the normalized maximum peak value of the waveform of the nth motor unit in M channels;

[0057] Step 3.3: Use Equation (7) to establish the regression relationship between ρ n and the twitch force peak value P n of the nth motor unit:

[0058] P n = a × ρ n + b (7)

[0059] In Equation (7), a and b are two regression coefficients, and their constraints are non - negative;

[0060] Step 3.4: Use Equation (8) to calculate the contraction time T n of the nth motor unit:

[0061]

[0062] In Equation (8), T L is the longest contraction time;

[0063] Step 3.5: Use Equation (9) to calculate the twitch force waveform k of the nth motor unit:

[0064]

[0065] In Equation (9), ω is the time - range vector;

[0066] Step 3.6: Execute Steps 3.3 to 3.5 for each of the N motor units respectively, so as to obtain the twitch force waveforms χ = [χ1, …, χ​n …, χ N T ;

[0067] For all motor units of the i-th muscle tendon unit, after convolving their firing sequences with the twitch force waveform respectively and then performing superposition processing, the original neural activation vector of the i-th muscle tendon unit is obtained After processing all k muscle tendon units, the original neural activation vectors of the k muscle tendon units are obtained

[0068] Step 3.7: Obtain the normalized neural activation vector of the i-th muscle tendon unit using Equation (10) Thus, k normalized neural activation vectors are obtained

[0069]

[0070] In Equation (10), is the maximum original activation value that the i-th muscle tendon unit can generate, that is, the peak value of the original neural activation vector obtained when setting ρ as a vector of all 1s;

[0071] Step 3.8: Process the residual electromyogram waveform matrix QResidual according to the non-negative matrix factorization in Step 2.2 to obtain the time-varying curve H Residual and the weight matrix Q Residual ;

[0072] Extract the activation positions of each column in the weight matrix Q Residual using the method in Step 2.5 where represents the i-th residual activation position;

[0073] Combine the residual activation position p Residual with the position of the electrode relative to the muscle tendon unit, and pair the n residual activation positions with the k muscle tendon unit positions one by one; the corresponding time-varying curve is passed through a low-pass filter and amplified to be used as the residual neural activation of the i-th muscle tendon unit Thus, k residual neural activations are obtained

[0074] Step 3.9: Obtain the neural activation u of the i-th muscle tendon unit using Equation (11) i , thus obtaining the neural activations u1, u2,..., u i ,…, u k ; ​

[0075]

[0076] In Equation (11), k and c are constants.

[0077] The said step 4.3 includes:

[0078] Step 4.3.1: Calculate the muscle tendon length vector at time t from the joint angle θ(t) at time t and the musculoskeletal geometric model and the moment arm matrix R(t) = [r1(t), r2(t), …, r i (t), …, r k (t)] T , where represents the muscle tendon length of the i-th muscle tendon unit at time t;

[0079] Step 4.3.2: Calculate the active contraction force F CE,i (t) of the i-th muscle tendon unit at time t, so as to obtain the active contraction forces of k muscle tendon units at time t;

[0080] Step 4.3.3: Calculate the passive contraction force F PE,i (t) of the i-th muscle tendon unit at time t, so as to obtain the passive contraction forces of k muscle tendon units at time t;

[0081] Step 4.3.4: Calculate the muscle force of the i-th muscle tendon unit at time t so as to obtain the muscle forces of k muscle tendon units at time t;

[0082] Step 4.3.4.1: Calculate the pennation angle φ i (t) of the i-th muscle tendon unit at time t using Equation (12):

[0083]

[0084] Step 4.3.4.2: Calculate the muscle force of the i-th muscle tendon unit at time t using Equation (13)

[0085]

[0086] The said step 4.3.2 includes:

[0087] Step 4.3.2.1: Calculate the muscle fiber length of the i-th muscle tendon unit at time t using Equation (14)

[0088]

[0089] Step 4.3.2.2: Calculate the normalized muscle fiber length of the $i$-th muscle tendon unit at time $t$ using Equation (15).

[0090]

[0091] In Equation (15), $\lambda$ is a constant.

[0092] Step 4.3.2.3: Calculate the length-related force of active contraction in the $i$-th muscle tendon unit at time $t$ using Equation (16).

[0093]

[0094] In Equation (16), $\sigma$ is a constant.

[0095] Step 4.3.2.4: Calculate the muscle fiber contraction velocity $v$ of the $i$-th muscle tendon unit at time $t$ using Equation (17). i (t):

[0096]

[0097] In Equation (17), $dt$ is the time resolution.

[0098] Step 4.3.2.5: Calculate the normalized muscle fiber contraction velocity of the $i$-th muscle tendon unit at time $t$ using Equation (18).

[0099]

[0100] Step 4.3.2.6: Calculate the velocity-related force of active contraction in the $i$-th muscle tendon unit at time $t$ using Equation (19).

[0101]

[0102] Step 4.3.2.7: Calculate the active contraction force $F$ of the $i$-th muscle tendon unit at time $t$ using Equation (20). CE,i (t):

[0103]

[0104] The said Step 4.3.3 includes:

[0105] Step 4.3.3.1: Calculate the passive normalized muscle fiber length of the $i$-th muscle tendon unit at time $t$ using Equation (21).

[0106]

[0107] Step 4.3.3.2: Calculate the length-related force of passive contraction in the $i$-th muscle tendon unit at time $t$ using Equation (22).

[0108]

[0109] Step 4.3.3.3: Calculate the passive contraction force $F_{(t)}$ of the $i$-th muscle tendon unit at time $t$ using Equation (23). PE,i (t):

[0110]

[0111] An electronic device according to the present invention includes a memory and a processor, characterized in that the memory is used to store a program for supporting the processor to execute the joint angle prediction method, and the processor is configured to execute the program stored in the memory.

[0112] A computer-readable storage medium according to the present invention, characterized in that a computer program is stored on the computer-readable storage medium, and the computer program executes the steps of the joint angle prediction method when run by a processor.

[0113] The present invention can scientifically and physiologically decode the motor unit activity, and achieve more accurate and robust joint angle prediction through a muscle-skeletal model that complies with physiological laws. The specific beneficial effects are as follows:

[0114] 1. In step two of the present invention, the two-dimensional waveform information of the motor unit is fully utilized, and the non-negative matrix factorization algorithm is innovatively combined to complete the robust automatic extraction of the activation position of the motor unit, realizing the effective positioning of the activation position of the motor unit. Then, an unsupervised clustering algorithm is used to automatically cluster the activation positions of the motor units into a specific number of categories, completing the classification of the motor unit activities, which can also be regarded as a method for tracking motor units.

[0115] 2. In step three of the present invention, the assignment of the motor unit to the specific muscle tendon unit is completed in combination with the physiological position, realizing the connection between the microscopic neural commands and the mesoscopic muscle tendon units. Then, a regression model and a twitch force model are combined to complete the decoding of the motor unit to the neural command, fully utilizing the decomposed motor unit information and considering the contributions between different motor units, realizing the scientific and physiological decoding of the underlying neural commands.

[0116] 3. In step three of the present invention, a unique neural activation normalization method is adopted to ensure the rationality and robustness of the estimated neural activation. In addition, a hybrid method of motor unit neural activation and residual neural activation is adopted in step three to ensure the integrity and rationality of the information used, providing a guarantee for achieving a higher-precision estimation effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0117] Figure 1 It is a schematic flowchart of the method provided by the embodiment of the present invention;

[0118] Figure 2 It is an experimental schematic diagram provided by the embodiment of the present invention;

[0119] Figure 3 It is a schematic diagram of motor unit localization provided by the embodiment of the present invention;

[0120] Figure 4 It is a schematic diagram of the joint angle estimation result provided by the embodiment of the present invention. Detailed implementation manners

[0121] In this embodiment, a joint angle prediction method based on motor units and neuromusculoskeletal models is as Figure 1 shown, and includes the following steps:

[0122] Step 1. Acquisition and preprocessing of electromyography data and motion data during joint movement:

[0123] Step 1.1. Use an electromyography measurement device and high-density electrodes with M channels to collect high-density surface electromyography signals at time t during joint movement at the muscle to be measured, denoted as x(t)=[x1(t), x2(t),…, x m (t),…, x M (t)] T , where x m (t) represents the electromyography data of the m-th channel at time t, so as to obtain a high-density surface electromyography signal X with a duration of T1; T represents transpose;

[0124] At the same time, use an action acquisition device to collect the angle information of S degrees of freedom in the joint, denoted as where represents the angle of the s-th degree of freedom in the joint at time t, so as to obtain the true joint angle with a duration of T1

[0125] Use an electromyography measurement device to collect the electromyography peak value when the subject performs a maximum voluntary contraction, denoted as MVC;

[0126] Step 1.2. Use a blind source separation algorithm to decompose the high-density surface electromyography signal X with a duration of T1 to obtain N motor unit firing sequences S=[s1, s2,…, s n ,…, s N T and its two-dimensional waveform W=[w 11 , w 12 ,…, w 1M ,…, w nm ​,…,w NM T and the residual electromyogram waveform matrix Residual, where s n represents the nth motor unit discharge sequence with a duration of T1, and w nm represents the waveform of the nth motor unit discharge sequence s n in the mth channel;

[0127] Specific implementations include:

[0128] (1) Recruit β subjects, guide the subjects to sit comfortably on a chair, place their arms naturally and comfortably on the chair handle, relax the shoulder joints, let the upper arms hang naturally, form a 0-degree angle with the body side, form a 90-degree angle at the elbows, and keep the forearms in a neutral position. Carefully clean the skin of the subjects' forearms with alcohol, and then attach the B-sheet array electrodes to the flexor and extensor muscle groups on the subjects' forearms. The radius of a single electrode contact in the array electrodes is radius, and the center-to-center spacing of the electrodes is interval. Use the electromyogram acquisition device with channels to collect the electromyogram signals of the flexor and extensor muscle groups of the subjects' forearms. The motion acquisition device can use an optical capture system, such as the Optitrack optical motion capture system. Specifically, place F marker points on the subjects' forearms, wrists, and hands to record the hand movements of the subjects, and collect the motion trajectories of the marker points through Optitrack, as Figure 2 shown. Exemplarily, it can be set that: β = 8, B = 2, radius = 2mm, interval = 10mm, F = 7.

[0129] (2) At the start of the experiment, synchronously collect the continuous electromyogram signals and the motion trajectories of the marker points when the subjects perform the wrist flexion and extension tasks, and use a trigger signal as a synchronization marker. The wrist flexion and extension tasks are divided into periodic motion and voluntary motion. In the periodic motion, the subjects are required to perform the wrist flexion-return to normal-extension-return to normal actions periodically, and the frequency of the whole process is controlled at about 1 / 8 Hz, and a total of five groups are completed; in the voluntary motion, no regulations are made on the motion direction, speed, and angle of the subjects' wrists, and the subjects can complete the wrist flexion and extension tasks in any direction, at any speed, and at any angle within about 20 s, as Figure 2 shown.

[0130] ​(3) Save the collected EMG and marker data offline. Scale the general wrist model in OpenSim software to match the body shape of a specific subject. Map the 3D motion of the markers to the wrist joint angles through an inverse kinematics tool, and upsample the angular data frequency to be consistent with the EMG frequency for further analysis. According to the trigger marker at the start of the experiment, further synchronize the EMG and angular data. Exemplarily, the collected high-density EMG signals can be decomposed using a stepwise variable peeling algorithm, and the obtained motor unit activities are saved offline for further analysis.

[0131] Step 2: Use the non-negative matrix factorization algorithm to extract the activation positions of motor units and cluster them:

[0132] Since the blind source separation algorithm requires sufficient spatial information, the electrodes are mostly high-density array electrodes with a large area. In practical applications, it is often unrealistic to place high-density electrodes of a specific shape on a specific muscle to only collect the activities of a target muscle. Therefore, high-density electrodes often collect the activities of multiple muscle tendon units, and the obtained motor unit activities are also from multiple muscle tendon units. To scientifically and physiologically decode motor unit information, it is necessary to locate and classify dozens of obtained motor units to determine their approximate activation positions. In the present invention, through the non-negative matrix factorization algorithm, the main activation positions are extracted from the two-dimensional waveform matrix of motor units, with good accuracy and robustness. At the same time, using an unsupervised clustering strategy, the activation positions of motor units are clustered to complete the classification of motor units, which conforms to the underlying physiological laws and successfully realizes the positioning and classification of motor units.

[0133] Step 2.1: Rearrange the two-dimensional waveform W of N M channels into N waveform matrices V1, V2, … V n , … V N , where V n represents the waveform matrix of the nth motor unit, each row of which represents a channel and each column represents a time sampling point. In this embodiment, the electrodes used are two square-shaped 64-channel array electrodes. Therefore, in Step 2.1, M = 64. In addition, the muscles that complete the wrist flexion and extension movements are simplified into 4 pieces, including two extensor muscles: extensor carpi radialis longus and extensor carpi ulnaris, and two flexor muscles: flexor carpi radialis and flexor carpi ulnaris. Since only one electrode is placed on the surface of each of the wrist flexor muscle group and extensor muscle group in this embodiment, the data of the two electrodes can be simply analyzed separately.

[0134] Step 2.2: Decompose the waveform matrix V n of the nth motor unit using the non-negative matrix factorization algorithm shown in Equation (1):

[0135]

[0136] In Equation (1), k is the number of muscle tendon units, Q n is the activation weight matrix of the nth motor unit, and H n is the time-varying curve matrix of the nth motor unit; in this embodiment, the k value of non-negative matrix factorization is set to 2, that is, it is considered that the data of one electrode comes from two muscle tendon units.

[0137] Step 2.3: For each row in H n , solve its peak value to obtain k time-varying curve peak values where δ i represents the peak value of the ith time-varying curve; for each column in Q n , solve the peak value to obtain k activation weight peak values Multiply the k time-varying curve peak values and the k activation weight peak values one by one in order to obtain k waveform peak values, and record the position index of the maximum value among the k waveform peak values as γ n .

[0138] Step 2.4: Take the γ n th column of the activation weight matrix Q n as the main activation pattern, and rearrange its M activation weights into the arrangement form of the high-density electrode.

[0139] Step 2.5: Set an activation weight threshold ε, obtain the two-dimensional positions greater than ε among the above M activation weights and take the average to obtain the activation position p n of the nth motor unit; in this embodiment, the activation weight threshold ε set in Step 2.5 is used to screen several channel positions with the largest weight values, and can be exemplarily set to 0.6 times the activation weight peak value. The specific extraction process is as Figure 3 shown.

[0140] Step 2.6: Obtain the activation positions p1,..., p n ..., p N of all N motor units according to the process from Step 2.2 to Step 2.5, and record them as P = [p1,..., p n ..., p N T , where p n represents the activation position of the nth motor unit firing sequence s n .

[0141] Step 2.7: Let the number of cluster centers be the same as the assumed number of muscle tendon units, denoted as k, and use the clustering algorithm to cluster the activation positions P of N motor units to obtain k cluster centers, denoted as C = [c1,..., c i ,..., c​k T and the class index g of each motor unit is g = [g1, …, g n …, g N T , where, c i represents the two-dimensional position of the i-th cluster center, and g n represents the class index of the n-th motor unit.

[0142] Step 3. Estimate neural activation through a regression model and a twitch force model:

[0143] After obtaining the activation positions of the motor units, it is necessary to further determine which muscle-tendon unit each motor unit belongs to, so as to determine the activities of each muscle-tendon unit. In the present invention, step 3.1 is used to determine the subordination of each motor unit; in addition, among the motor unit information obtained through surface electromyogram decomposition, there are motor unit firing information and waveform information, both of which contain rich spatio-temporal activity information of the motor units. Therefore, in the process of decoding the neural activation curve of the muscle-tendon unit from the motor unit activities, it is necessary to make full use of the firing information and waveform information of each motor unit, and more importantly, to consider the differences between different motor units; the present invention completes the calculation of the neural activation curve from the motor unit activities to the muscle-tendon unit through steps 3.2 to 3.9, improving the scientificity and physiological nature of motor unit decoding.

[0144] Step 3.1. Pair the k categories with the k muscle-tendon units one by one according to the cluster center position C and the position of the electrode relative to the muscle-tendon unit; according to the class index g of the motor unit, assign the N motor unit firing sequences S to the k muscle-tendon units one by one; in this embodiment, the subordination of each type of motor unit can be determined according to the activation position of each type of motor unit on the electrode and the position of the muscle-tendon unit relative to the electrode. For example, the anatomical position of the flexor carpi radialis is on the radial side, and relative to the position where the electrode is placed, it is at the position of the 32nd channel from the back in the electrode array. Therefore, the type of motor unit with the activation position at the 32nd channel from the back belongs to the motor unit activity of the flexor carpi radialis.

[0145] Step 3.2. Calculate the maximum peak values of the waveforms of the N motor units in the M channels from W respectively, so as to obtain N maximum peak values, and normalize these N maximum peak values using MVC. The N normalized maximum peak values are denoted as ρ = [ρ1, ρ2, …, ρ n , …, ρ N T , ρ n ​​​represents the normalized maximum peak value of the waveform of the nth motor unit in M channels; in this embodiment, it is necessary to normalize the maximum value of the waveform with the electromyogram value MVC during maximum voluntary contraction. This value needs to require the subject to perform three maximum voluntary contractions at the beginning of the experiment and record their high-density electromyogram signals, and then extract the peak value of the electromyogram during the three contractions as the electromyogram value of the maximum voluntary contraction.

[0146] Step 3.3. Establish ρ using Equation (2) n and the twitch force peak value P of the nth motor unit n between the regression relationship:

[0147] P n = a × ρ n + b (2)

[0148] In Equation (2), a and b are two regression coefficients, and their constraints are non-negative; in this embodiment, the linear regression model establishes the relationship between the peak value ρ of the motor unit action potential waveform and the twitch force amplitude P value, that is, the larger the peak value ρ of the waveform, the larger its P value. This actually establishes a connection between the size of the motor unit and the amplitude of the waveform. Under the current decomposition conditions, since it is very difficult to obtain the depth information of the motor unit, we can simply think that the larger the waveform of the motor unit, the larger its twitch amplitude and the greater its contribution to movement. This is a simplified idea. In previous model-free studies, it has been proven that this simplification can bring positive effects and enable the model to obtain more physiological information.

[0149] Step 3.4. Calculate the contraction time T of the nth motor unit using Equation (3) n :

[0150]

[0151] In Equation (3), T L is the longest contraction time.

[0152] Step 3.5. Calculate the twitch force waveform χ of the nth motor unit using Equation (4) n :

[0153]

[0154] In Equation (4), ω is the time range vector; in this embodiment, the twitch force model adopted has been verified by multiple previous studies and has good reliability and physiological interpretability. Among them, the longest contraction time can be set T L = 90ms.

[0155] Step 3.6: Perform the above Steps 3.3 to 3.5 on each of the N motor units, thereby obtaining the twitch force waveforms χ = [χ1, …, χ n …, χ N T ; For all the motor units of the i-th muscle-tendon unit, convolve their firing sequences with the twitch force waveforms respectively, and then perform superposition processing, thereby obtaining the original neural activation vector of the i-th muscle-tendon unit Perform the above process on all k muscle-tendon units, thereby obtaining the original neural activation vectors of the k muscle-tendon units

[0156] Step 3.7: Use Equation (5) to obtain the normalized neural activation vector of the i-th muscle-tendon unit Thereby obtaining k normalized neural activation vectors

[0157]

[0158] In Equation (5), is the maximum original activation value that the i-th muscle-tendon unit can generate, that is, the peak value of the original neural activation vector obtained when ρ is set as a vector of all 1s; since the neural activation curve is strictly limited between 0 and 1, a new normalization method is proposed in Step 3.7. This method divides the calculated neural curve by the maximum neural activation value that this muscle-tendon unit can generate, that is, the normalized ρ value of all motor units is 1. The reason for adopting this method is that manually defining the maximum neural activation value is too subjective, and other normalization methods such as min-max normalization are not applicable in this case. Therefore, the normalization method proposed by the present invention is the result of comprehensively considering rationality and accuracy.

[0159] Step 3.8: Process the residual EMG waveform matrix QResidual according to the non-negative matrix factorization in Step 2.2 to obtain the time-varying curve H Residual and the weight matrix Q Residual ;

[0160] Use the method in Step 2.5 to extract the activation positions of each column in the weight matrix Q Residual wherein represents the i-th residual activation position; Combined with the residual activation position p

[0161] and the position of the electrode relative to the muscle-tendon unit, pair the n residual activation positions with the k muscle-tendon unit positions one by one; the one paired with the i-th muscle-tendon unit Residual ​​After passing through a low-pass filter and being amplified, the corresponding time-varying curve serves as the residual neural activation of the i-th muscle-tendon unit. Thus, k residual neural activations are obtained. In this embodiment, the cut-off frequency of the low-pass filter in step 3.8 can be set to 4 Hz.

[0162] Step 3.9: Use Equation (6) to obtain the neural activation u of the i-th muscle-tendon unit. i , thus obtaining the neural activations u1, u2, …, u of k muscle-tendon units. i , …, u k ;

[0163]

[0164] In Equation (6), κ and c are constants.

[0165] In this embodiment, since the motor unit information is not available at all times during the angular task. For example, the number of motor units at small angles is very small or even non-existent. Therefore, it is necessary to superimpose the residual EMG information. So, in the case where there is no firing or the number of motor unit firings is extremely small, we use the EMG residual to calculate the neural activation curve. And because their processing methods are different and there is a difference in amplitude, it is necessary to use a scaling factor to unify them; κ in Equation (6) is an empirical value determined through preliminary experiments and is set to 3 in this embodiment.

[0166] Step 4: In this embodiment, based on the phenomenological model - Hill model, a musculoskeletal model for wrist flexion and extension movements is modeled to complete the solution from neural activation to joint angle. This model includes the following modules: muscle activation dynamics, muscle-tendon model, joint geometry model, and dynamics model, which can gradually transform the extracted neural activation curve into the macroscopic limb movement angle.

[0167] Step 4.1: Nonlinearize U using Equation (7) to obtain the corresponding muscle activation A:

[0168]

[0169] In Equation (7), z is a nonlinear factor, A = [a1, a2, …, a i , …, a k T , where, a i = [a i (1), a i (2), …, a i (t), …, a i (T1)], a i ​(t) represents the muscle activation of the i-th muscle-tendon unit at time t; in this embodiment, the non-linear factor z needs to be set with a range constraint, and can be set to -3 to 0 exemplarily.

[0170] Step 4.2, initialize t = 0;

[0171] Set the estimated joint angle θ(t) = [θ1(t), θ2(t), …, θ s (t), …, θ S (t)] at time t, where θ s (t) represents the angle of the s-th degree of freedom in the joint estimated at time t;

[0172] Set the muscle fiber contraction velocity v(t) = [v1(t), v2(t), …, v i (t), …, v k (t)] T as the zero vector, where v i (t) represents the muscle fiber contraction velocity of the i-th muscle-tendon unit at time t;

[0173] Given the parameter vector h of the musculoskeletal model:

[0174]

[0175] In Equation (8), φ o,i , and are respectively the maximum isometric force, optimal muscle fiber length, optimal pennation angle, tendon length and length scaling factor of the i-th muscle-tendon unit.

[0176] Step 4.3, calculate the muscle forces of k muscle-tendon units at time t using the Hill muscle model where, represents the muscle force of the i-th muscle-tendon unit at time t: A series of modeling formulas in Step 4.3 can convert the muscle activation curve into muscle force, and these modeling formulas are generally recognized internationally and describe the phenomenological response of muscles when stimulated to contract.

[0177] Step 4.3.1, calculate the muscle-tendon length vector at time t from the joint angle θ(t) at time t and the musculoskeletal geometric model and the moment arm matrix R(t) = [r1(t), r2(t), …, r i (t), …, r k (t)] T , where, Represents the muscle-tendon length of the i-th muscle-tendon unit at time t; the musculoskeletal geometric model in this embodiment can be modeled using OpenSim software. OpenSim is an open source software developed by Stanford University, and there are many public resources for researchers to use. For example, the musculoskeletal geometric model can be selected as the general wrist model that comes with OpenSim.

[0178] Step 4.3.2: Calculate the active contraction force F of the ith muscle-tendon unit at time t CE,i (t), thus obtaining the active contraction force of k muscle-tendon units at time t:

[0179] Step 4.3.2.1. Calculate the muscle fiber length of the i-th muscle-tendon unit at time t using formula (9):

[0180]

[0181] Step 4.3.2.2: Calculate the normalized muscle fiber length of the ith muscle-tendon unit at time t using formula (10):

[0182]

[0183] In formula (10), λ is a constant; in this embodiment, λ can be set to 0.15.

[0184] Step 4.3.2.3. Calculate the length-dependent force of active contraction in the i-th muscle-tendon unit at time t using equation (11):

[0185]

[0186] In formula (11), σ is a constant; in this embodiment, σ can be set to 0.45.

[0187] Step 4.3.2.4: Calculate the muscle fiber contraction velocity v of the ith muscle-tendon unit at time t using equation (12): i (t):

[0188]

[0189] In formula (12), dt is the time resolution;

[0190] Step 4.3.2.5: Calculate the normalized muscle fiber contraction velocity of the ith muscle-tendon unit at time t using formula (13):

[0191]

[0192] Step 4.3.2.6: Calculate the velocity-related force of active contraction in the \(i\)-th muscle-tendon unit at time \(t\) using Equation (14).

[0193]

[0194] Step 4.3.2.7: Calculate the active contraction force \(F_{i}(t)\) in the \(i\)-th muscle-tendon unit at time \(t\) using Equation (15). CE,i (t):

[0195]

[0196] Step 4.3.3: Calculate the passive contraction force \(F_{i}^{p}(t)\) of the \(i\)-th muscle-tendon unit at time \(t\), so as to obtain the passive contraction forces of \(k\) muscle-tendon units at time \(t\): PE,i (t):

[0197] Step 4.3.3.1: Calculate the passive normalized muscle fiber length of the \(i\)-th muscle-tendon unit at time \(t\) using Equation (16).

[0198]

[0199] Step 4.3.3.2: Calculate the length-related force of passive contraction in the \(i\)-th muscle-tendon unit at time \(t\) using Equation (17).

[0200]

[0201] Step 4.3.3.3: Calculate the passive contraction force \(F_{i}^{p}(t)\) of the \(i\)-th muscle-tendon unit at time \(t\) using Equation (18). PE,i (t):

[0202]

[0203] Step 4.3.4: Calculate the total force of the \(i\)-th muscle-tendon unit at time \(t\).

[0204] Step 4.3.4.1: Calculate the pennation angle \(\varphi_{i}(t)\) of the \(i\)-th muscle-tendon unit at time \(t\) using Equation (19). i (t):

[0205]

[0206] Step 4.3.4.2: Calculate the muscle force of the \(i\)-th muscle-tendon unit at time \(t\) using Equation (20).

[0207]

[0208] Step 4.4: Obtain the torques τ(t) of S degrees of freedom in the joint at time t using Equation (21):

[0209]

[0210] In Equation (21), τ(t) = [τ1(t), τ2(t), …, τ s (t), …, τ S (t)], where τ s (t) represents the torque of the s-th degree of freedom in the joint at time t; r i (t) represents the moment arm vector of the i-th muscle-tendon unit at time t, and r i (t) = [r i1 (t), r i2 (t), …, r is (t), …, r iS (t)], and r is (t) is the moment arm of the i-th muscle-tendon unit at the s-th degree of freedom in the joint.

[0211] Step 4.5: Solve the state equation shown in Equation (22) using numerical integration to obtain the estimated joint angle θ(t + 1) at time t + 1;

[0212]

[0213] In Equation (22), M is the inertia matrix, C is the damping matrix, K is the stiffness matrix, and G is the gravity matrix; represents the first derivative of the estimated joint angle at time t, represents the second derivative of the estimated joint angle at time t;

[0214] In this embodiment, the inertia, damping, stiffness, and gravity terms need to be calculated according to the musculoskeletal model. Exemplarily, in this experiment, M = mι 2 , C = 0.3, K = 1.5088, G = 0, where m is the mass of the hand and ι is the length from the wrist rotation center to the center of mass of the hand.

[0215] In this embodiment, the numerical integration method used to solve Equation (22) is the first-order Euler method. Similarly, other numerical integration methods such as the fourth-order Runge-Kutta method can be adopted.

[0216] Step 4.6: After assigning t + 1 to t, if t ≠ T1, return to Step 4.3 and execute sequentially to solve the estimated joint angle at the next moment;

[0217] Step 5: Optimize the extended parameter vector in the musculoskeletal model based on the global heuristic search algorithm:

[0218] Step 5.1. Initialize the extended parameter vector

[0219]

[0220] Step 5.2. Use the extended parameter vector Process the high-density surface electromyogram signal X with a duration of T1 according to the process of Steps 2 - 4 to obtain the estimated joint angles with a duration of T1, denoted as Θ = [θ(0), θ(1), …, θ(t), …, θ(T1)];

[0221] Step 5.3. Using the true joint angles as a reference, within the set constraints, use a global optimization algorithm to optimize and solve the extended parameter vector so that the root mean square error between the estimated joint angle Θ and the true joint angle is minimized, thereby obtaining the optimized extended parameter vector Use the determined neuromusculoskeletal model to achieve the prediction of joint angles.

[0222] In this embodiment, the global optimization algorithm adopted is the genetic algorithm. Similarly, other global optimization algorithms such as the simulated annealing algorithm and the particle swarm algorithm can be adopted; the parameters to be optimized in the global optimization algorithm need to set initial values and ranges. In this embodiment, the initial parameters of all optimized muscle tendon units are the same as those of the general musculoskeletal model in the OpenSim software, and the optimization range can set a certain interval within the physiological range. Exemplarily, φ o,i can be respectively constrained within ±50%, ±5%, ±5%, ±5% of the initial value, is constrained between 0.9 and 1.2, and c is constrained between 0 and 3;

[0223] After the parameter optimization is completed, the optimized extended parameter vector can be obtained. Apply the parameter values in to the neuromusculoskeletal model in Steps 2 to 4, thereby obtaining the subject-specific neuromusculoskeletal model. In the joint angle prediction stage, only by collecting the high-density surface electromyogram signal of the subject, the above neuromusculoskeletal model can be used for accurate joint angle prediction, as Figure 4 shown.

[0224] In summary, the present invention can effectively decode the underlying neural instruction information into macroscopic joint movements, and the intermediate process is transparently visible, which is reasonable and physiologically meaningful. It has better robustness and physiological interpretability compared with machine learning or deep learning methods, and has more detailed and accurate advantages compared with the method using macroscopic electromyogram envelopes. The proposed motor unit decoding framework of the present invention can complete the localization of the activation position of the motor unit and the allocation to specific muscle tendon units, and complete the functional classification and fuzzy tracking of the motor unit in a scientific and physiological manner. In addition, the present invention makes full use of the neural instruction information obtained by decomposition, and both the motor unit firing sequence and waveform information are used to estimate the neural activation curve, which can more comprehensively explore the underlying neural instructions. The estimation method of the neural activation curve is more novel, detailed and effective compared with traditional methods. The present invention proposes a new method for predicting joint angles based on motor units and neuromusculoskeletal models, providing a new solution and idea for the difficulties in motor unit function decoding and the unclear representation of the relationship between microscopic neural instructions and macroscopic limb movements faced in current practical motion control research.

[0225] In this embodiment, an electronic device includes a memory and a processor. The memory is used to store a program that supports the processor to execute the above method, and the processor is configured to execute the program stored in the memory.

[0226] In this embodiment, a computer-readable storage medium stores a computer program, and when the computer program is run by a processor, it executes the steps of the above method.

Claims

1. A joint angle prediction method based on motor units and neuromusculoskeletal models, characterized in that, It is carried out according to the following steps: Step 1: Acquisition and preprocessing of EMG data and motion data during joint movement: Step 1.1: Use an electromyogram measurement device and high-density electrodes with M channels to collect high-density surface electromyogram signals at time t during joint movement at the muscle to be measured, denoted as x(t) = [x1(t), x2(t),..., x m (t),..., x M (t)] T , where x m (t) represents the electromyogram data of the m-th channel at time t, thereby obtaining a high-density surface electromyogram signal X with a duration of T1; T represents the transpose; Meanwhile, an action acquisition device is used to collect the angular information of S degrees of freedom in the joint, denoted as where represents the angle of the s-th degree of freedom in the joint at time t, thereby obtaining the true joint angle with a duration of T1 Use an EMG measurement device to collect the peak EMG value when the subject performs a maximum voluntary contraction, denoted as MVC; Step 1.

2. Decompose the high-density surface electromyogram signal X with a duration of T1 by using a blind source separation algorithm to obtain N motor unit discharge sequences S = [s1, s2,..., s n ,..., s N T and its two-dimensional waveform W = [w 11 , w 12 ,..., w 1M ,..., w nm ,..., w NM R and the residual electromyogram waveform matrix Residual, where s n represents the nth motor unit discharge sequence with a duration of T1, and w nm represents the waveform of the nth motor unit discharge sequence s n in the mth channel;​​ Step 2: Extract the motor unit activation positions and cluster them using the non-negative matrix factorization algorithm: Extract the activation positions of N motor units from the two-dimensional waveform W of N M channels, denoted as P = [p1,..., p n ... , P N T , where p n represents the activation position of the firing sequence s n of the nth motor unit;​ Let the number of cluster centers be the same as the assumed number of muscle tendon units, denoted as k. Then, use the clustering algorithm to cluster the activation positions P of N motor units to obtain k cluster centers denoted as C = [c1,..., c i ,...c k T and the class indices g of each motor unit, g = [g1,..., g n ..., g N T , where c i represents the two-dimensional position of the i-th cluster center, and g n represents the class index of the n-th motor unit;​​ Step 3: Estimate neural activation through a regression model and a twitch force model: Estimate the neural activations of k muscle-tendon units from S and W, denoted as U = [u1, u2,..., u i ,..., u k T , where u i represents the neural activation of the i-th muscle-tendon unit with a duration of T1, and u i = [u i (1), u i (2),..., u i Lt),..., u i (T1)], where u i (t) represents the neural activation of the i-th muscle-tendon unit at time t;​ Step 4: Complete the calculation from neural activation to joint angle based on a musculoskeletal model: Step 4.1: Nonlinearize U using Equation (1) to obtain the corresponding muscle activation A: In formula (1), z is a non-linear factor, A = [a1, a2,..., a i ,..., a k R , where a i = [a i (1), a i (2),..., a i (t),..., a i (T1)], a i (t) represents the muscle activation of the i-th muscle tendon unit at time t;​ Step 4.2: Initialize t = 0; Set the estimated joint angle θ(t) at the t-th moment as θ(t) = [θ1(t), θ2(t),..., θ s (t),..., θ S (t)], where θ s (t) represents the angle of the s-th degree of freedom in the estimated joint at the t-th moment; Set the muscle fiber contraction velocity v(t) at the t-th moment = [v1(t), v2(t),..., v i (t),..., v k (t)] T to be the zero vector, where v i (t) represents the muscle fiber contraction velocity of the i-th muscle tendon unit at the t-th moment; Given the parameter vector h of the musculoskeletal model: In Equation (2), φ o，i , and are respectively the maximum isometric force, optimal muscle fiber length, optimal pennation angle, tendon length, and length scaling factor of the i-th muscle tendon unit; Step 4.3: Calculate the muscle forces of k muscle-tendon units at time t using the Hill muscle model where represents the muscle force of the i-th muscle-tendon unit at time t; Step 4.4: Use Equation (3) to obtain the torque τ(t) of S degrees of freedom in the joint at time t: In Equation (3), τ(t) = [τ1(t), τ2(t),..., τ s (t),..., τ S (t)], where τ s (t) represents the torque of the s-th degree of freedom in the joint at time t; r i (t) represents the moment arm vector of the i-th muscle-tendon unit at time t, and r i (t) = [r i1 (t), r i2 (t),..., r is (t),..., r iS (t)], where r is (t) is the moment arm of the i-th muscle-tendon unit at the s-th degree of freedom in the joint; Step 4.5: Solve the state equation shown in Equation (4) using numerical integration to obtain the estimated joint angle θ(t + 1) at time t + 1; In Equation (4), M is the inertia matrix, C is the damping matrix, K is the stiffness matrix, and G is the gravity matrix; represents the first-order differential of the estimated joint angle at time t, represents the second-order differential of the estimated joint angle at time t; Step 4.6: After assigning t + 1 to t, if t ≠ T1, return to Step 4.3 and execute sequentially to solve the estimated joint angle at the next moment; Step 5: Optimize the extended parameter vector in the musculoskeletal model based on a global heuristic search algorithm: Step 5.1, initialize the extended parameter vector In Equation (5), a and b represent two regression coefficients, and c is an amplitude scaling factor; Step 5.

2. Using the extended parameter vector Process the high-density surface electromyogram signal X with a duration of T1 according to the process of Steps 2 - 4 to obtain the estimated joint angles with a duration of T1, denoted as Θ = [θ(0), θ(1),..., θ(t),..., θ(T1)]; Step 5.

3. Using the true joint angle as a reference, within the set constraints, use a global optimization algorithm to optimize and solve the extended parameter vector so that the root mean square error between the estimated joint angle Θ and the true joint angle is minimized, thereby obtaining the optimized extended parameter vector Use the neuromusculoskeletal model determined thereby to predict the joint angle.

2. The joint angle prediction method based on a motor unit and a neuromusculoskeletal model according to claim 1, wherein The extraction and clustering of the motor unit activation positions in Step 2 are carried out according to the following steps: Step 2.1: Rearrange the N two-dimensional waveforms W with M channels into N waveform matrices V1, V2,... V n ,... V N , where V n represents the waveform matrix of the nth motor unit, each row of which represents a channel and each column represents a time sampling point; Step 2.2: Decompose the waveform matrix V of the nth motor unit by using the non-negative matrix factorization algorithm shown in Equation (6) n as follows: In Equation (6), k is the number of muscle tendon units, Q n is the activation weight matrix of the nth motor unit, and H n is the time-varying curve matrix of the nth motor unit; Step 2.

3. For each row in H n , solve for its peak value to obtain k peak values of time-varying curves where δ i represents the peak value of the i-th time-varying curve; For Q n For each column in it, solve for its peak value to obtain k activation weight peak values Multiply the k time-varying curve peaks and the k activation weight peaks one by one in order to obtain k waveform peaks, and record the position index of the maximum value among the k waveform peaks as γ n ; Step 2.4: Take the γ-th column of the activation weight matrix Q n as the main activation pattern, and rearrange its M activation weights into the arrangement form of the high-density electrodes; n ​ Step 2.5: Set an activation weight threshold ε, obtain the two-dimensional positions among the M activation weights that are greater than ε, and take the average to obtain the activation position p of the nth motor unit n ; Step 2.

6. Obtain the activation positions p1,..., p of all N motor units according to the process from Step 2.2 to Step 2.5 n ...,p N ; Step 2.7: Use a clustering algorithm to cluster the activation positions P of N motor units to obtain k cluster centers denoted as C and the class indices g of each motor unit.

3. The joint angle prediction method based on motor units and neuromusculoskeletal models according to claim 1, characterized in that Step 3 includes: Step 3.1: Pair the k classes with k muscle-tendon units one by one according to the cluster center position C and the position of the electrode relative to the muscle-tendon unit; according to the class index g of the motor unit, assign the firing sequences S of N motor units to the k muscle-tendon units one by one; Step 3.2: Calculate the maximum peak values of the waveforms of N motor units in M channels from W respectively, so as to obtain N maximum peak values, and normalize the N maximum peak values using MVC to obtain N normalized maximum peak values, denoted as ρ = [ρ1, ρ2,... n ,... N , ρ T , where ρ n represents the normalized maximum peak value of the waveform of the nth motor unit in M channels; Step 3.

3. Establish the regression relationship between ρ n and the peak twitch force P n of the nth motor unit: P n = a × ρ n + b (7) In Equation (7), a and b are two regression coefficients, and their constraints are non-negative; Step 3.

4. Calculate the contraction time T of the nth motor unit using Equation (8). n :[[]]END]] In formula (8), T L is the longest contraction time; Step 3.

5. Calculate the twitch force waveform χ of the nth motor unit using Equation (9). n : In Equation (9), ω is a time range vector; Step 3.6: Perform Steps 3.3 to 3.5 on each of the N motor units, thereby obtaining the twitch force waveforms χ = [χ1,... χ n ... χ N T ;​ For all motor units of the i-th muscle tendon unit, their firing sequences are respectively convolved with the twitch force waveforms and then superimposed to obtain the original neural activation vector of the i-th muscle tendon unit. After processing all k muscle tendon units, the original neural activation vectors of the k muscle tendon units are obtained. Step 3.7: Obtain the normalized neural activation vector of the \(i\)-th muscle tendon unit using Equation (10). Thus, \(k\) normalized neural activation vectors are obtained. In Equation (10), is the maximum original activation value that the i-th muscle tendon unit can generate, that is, the peak value of the original neural activation vector obtained when ρ is set as a vector of all 1s; Step 3.8: Process the residual EMG waveform matrix Residual according to the non-negative matrix factorization in Step 2.2 to obtain the time-varying curve H Residual and the weight matrix Q Residual ; Extract the activation positions of each column in the weight matrix Q using the method in step 2.5 Residual wherein represents the i-th residual activation position ​ Combined with the residual activation position p Residual With the position of the electrode relative to the muscle-tendon unit, pair the k residual activation positions with the k muscle-tendon unit positions one by one; the corresponding time-varying curve is passed through a low-pass filter and amplified, and used as the residual neural activation of the i-th muscle-tendon unit Thus, k residual neural activations are obtained Step 3.9: Obtain the neural activation \(u\) of the \(i\)-th muscle-tendon unit using Equation (11) i , thereby obtaining the neural activations \(u_1, u_2,\cdots, u\) i , \(\cdots, u\) k ; In Equation (11), κ and c are constants.

4. The joint angle prediction method based on a motor unit and a neuromusculoskeletal model according to claim 1, wherein Step 4.3 includes: Step 4.3.1: Calculate the muscle tendon length vector at time t from the joint angle θ(t) at the t-th moment and the musculoskeletal geometric model and the moment arm matrix R(t) = [r1(t), r2(t),..., r i (t),..., r k (t)] T , where represents the muscle tendon length of the i-th muscle tendon unit at time t; Step 4.3.2, calculate the active contraction force F CE，i (t) of the i-th muscle-tendon unit at time t, so as to obtain the active contraction forces of k muscle-tendon units at time t; Step 4.3.3: Calculate the passive contraction force F PE，i (t) of the i-th muscle tendon unit at time t, so as to obtain the passive contraction forces of k muscle tendon units at time t; Step 4.3.4, calculate the muscle force of the i-th muscle tendon unit at time t Thus, the muscle forces of k muscle tendon units at time t are obtained; Step 4.3.4.

1. Calculate the pennation angle φ i (t) of the i-th muscle-tendon unit at time t using Equation (12): Step 4.3.4.2: Calculate the muscle force of the $i$-th muscle-tendon unit at time $t$ using Equation (13).

5. The calculation of the muscle force of the i-th muscle-tendon unit at time t according to claim 4, characterized in that, Step 4.3.2 includes: Step 4.3.2.

1. Calculate the muscle fiber length of the $i$-th muscle tendon unit at time $t$ using Equation (14). Step 4.3.2.

2. Calculate the normalized muscle fiber length of the $i$-th muscle tendon unit at time $t$ using Equation (15). In Equation (15), λ is a constant; Step 4.3.2.3: Calculate the length-related force of active contraction in the $i$-th muscle-tendon unit at time $t$ using Equation (16). In Equation (16), σ is a constant; Step 4.3.2.4, calculate the muscle fiber contraction velocity v i (t) of the i-th muscle tendon unit at time t using Equation (17): In Equation (17), dt is the time resolution; Step 4.3.2.

5. Calculate the normalized muscle fiber contraction velocity of the $i$-th muscle tendon unit at time $t$ using Equation (18). Step 4.3.2.6, calculate the velocity-related force of active contraction in the $i$-th muscle-tendon unit at time $t$ using Equation (19) Step 4.3.2.7, calculate the active contractile force F CE，i (t) in the i-th muscle-tendon unit at time t using Equation (20):

6. The calculation of the muscle force of the i-th muscle-tendon unit at time t according to claim 4, characterized in that Step 4.3.3 includes: Step 4.3.3.1: Calculate the passive normalized muscle fiber length of the $i$-th muscle tendon unit at time $t$ using Equation (21). Step 4.3.3.

2. Calculate the length-related force of passive contraction in the i-th muscle tendon unit at time t using Equation (22). Step 4.3.3.

3. Calculate the passive contractile force F of the i-th muscle-tendon unit at time t using Equation (23): PE，i (t):

7. An electronic device, comprising a memory and a processor, characterized in that, The memory is used to store a program that supports the processor to execute the joint angle prediction method described in any one of claims 1 - 6, and the processor is configured to execute the program stored in the memory.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is run by the processor, it executes the steps of the joint angle prediction method described in any one of claims 1 - 6.

Images