A method for upper limb rehabilitation assessment of stroke patients based on electroencephalomyogram signals

Abstract

The present invention discloses a method for upper limb rehabilitation assessment of stroke patients based on electroencephalomyogram signals, which specifically includes the following steps: signal acquisition and preprocessing; electroencephalomyogram signal decomposition; multi-scale processing; transfer spectral entropy calculation; upper limb rehabilitation level assessment. Based on the Spearman correlation coefficient, the fitting relationship between the significant area value of MSTSE of each patient and the Brunnstrom assessment scale is mapped, the FCMC characteristic indexes and change trends that are significantly correlated with the levels of the Brunnstrom assessment scale are analyzed, and the quantitative assessment of the disease levels of stroke patients is carried out. The technical solution of the present invention uses the binary variational mode decomposition-multi-scale transfer spectral entropy method to characterize the FCMC coupling relationship between the motor cortex of the brain and the corresponding muscle tissues, and combines with the Brunnstrom assessment scale to realize the objective quantitative assessment of the upper limb rehabilitation levels of stroke patients.

Description

Technical Field

[0001] The present invention relates to the technical field of bioelectrical signal processing, and particularly relates to an upper limb rehabilitation assessment method for stroke patients based on electroencephalogram and electromyogram signals. Background Art

[0002] In recent years, the incidence risk of stroke patients ranks first globally. Each year, more than 20 million patients suffer from varying degrees of limb dysfunction caused by cerebrovascular injuries, which has a serious negative impact on people's normal life, work, family and society. Effective targeted rehabilitation training is beneficial for patients to recover related motor functions, but accurate post-stroke rehabilitation assessment is a prerequisite for various rehabilitation trainings. Currently, traditional clinical assessment methods are mostly based on clinical assessment scales by doctors, and the level of assessment is completed according to the degree of action execution of patients. This method has the following disadvantages: the assessment results are greatly affected by the subjective factors of doctors. Different doctors may give different judgment results for the same patient based on different assessment scales and their experiences; some subtle manifestations of patients are not timely discovered by doctors, resulting in inaccurate assessment results.

[0003] With the development of medical technology, as well as computer and electronic technologies, clinical diagnosis and rehabilitation methods after stroke have gradually become diverse. Electroencephalogram (EEG) and surface electromyogram (sEMG) can quickly characterize the changes in the functions of the damaged brains and muscles of patients, and are widely used in clinical practice due to their low cost. When the human body autonomously executes an action, the EEG formed by the cerebral cortex first reaches the motor unit through the central nervous system. The motor unit controls muscle fibers to complete muscle contraction, generating action potentials on the muscle fibers to achieve action execution. This process reflects the information on the control of limb movement by the brain and the response information of muscles to the control of movement by the brain, which is called functional cortico-muscular coupling (FCMC). A large number of studies have shown that synchronizing EEG and sEMG information can characterize the functional coupling between the brain and its innervated muscles, reflecting the functional state of the motor nervous system. Stroke causes damage to the cerebral cortex and nerve pathways, resulting in local neurological deficits in the brain and the inability to send accurate and complete nerve impulses to the peripheral nervous system, thereby causing motor dysfunction.

[0004] There are currently various methods for studying FCMC, such as Granger Causality (GC), Transfer Entropy (TE), generalized partial directed coherence (gPDC), etc. GC can effectively calculate the causal relationship between signal sequences and can reflect the directionality between cortical and muscle coupling. However, GC can only measure the linear relationship between signals, and both EEG and sEMG are random signals. Therefore, the non-linear relationship between signals cannot be well characterized. The TE algorithm quantitatively analyzes the non-linear directed coupling relationship between signals and is widely used in FCMC research. However, TE can only capture the information flow relationship in one direction, and for bidirectional or multi-directional information flow in complex systems, TE lacks an effective explanation. The gPDC algorithm can distinguish the directional connection between signals, determine the direction of information transfer, and is applicable to the causal relationship analysis of multiple time series in complex systems. However, this algorithm is relatively sensitive to noise and is greatly affected by parameters.

[0005] Since both electroencephalogram and electromyogram signals are random signals and contain complex features at different scales, the use of decomposition algorithms is relatively common in FCMC, such as: wavelet packet, empirical mode decomposition (EMD), variational mode decomposition (VMD), etc. The wavelet packet algorithm is often affected by the selection of basis functions, and improper selection will affect the decomposition effect and analysis results. The EMD algorithm has an end effect, is inherently uncertain, and lacks mathematical theory support. The VMD algorithm has high resolution and strong self-adaptability, but it is greatly affected by initial parameters and can only process single-channel signals. Summary of the Invention

[0006] The technical solution of the present invention to solve the above technical problems is to provide an upper limb rehabilitation assessment method for stroke patients based on electroencephalogram and electromyogram signals, including the following steps:

[0007] Signal acquisition and preprocessing;

[0008] Decomposition of electroencephalogram and electromyogram signals:

[0009] Construct a complex signal, combine the electroencephalogram signal x(t) and the electromyogram signal y(t) into a complex form, defined as: z(t) = x(t) + jy(t), where j is the imaginary unit;

[0010] Direction selection, select a set of directions to project z(t), and select multiple different directions θ d , where θ dis a uniformly distributed angle on the unit circle. In the case of selecting D directions, each direction θ d can be expressed as: d = 0, 1, 2, …, D - 1;

[0011] Projection calculation: Project the complex signal z(t) onto a certain direction θ, and calculate the projection p θ (t) of z(t) in this direction. p θ (t) = Re(z(t)e -jθ ), where Re represents taking the real part, and θ is the angle of the projection direction;

[0012] Parameter initialization: Initialize the parameters required for VMD, including the number of modes K, the center frequencies w k of each mode, the Lagrange multiplier λ, the mode function u k , the penalty parameter α, and the convergence tolerance tol;

[0013] Define the objective function, requiring that each mode u k (t) has a finite bandwidth and center frequency. The constraint condition is that the sum of all modes is equal to the input signal p θ (t), and the sum of the estimated bandwidths of each mode is minimized. The constrained variational model is as follows:

[0014]

[0015] s.t. ∑ k u k = p θ (t);

[0016] In the formula, {w k} represents the frequency centers of each intrinsic mode function (IMFs), {u k} are the K finite-bandwidth IMFs components obtained after decomposition, is for taking partial derivatives, * is the convolution operator, δ(t) is the Dirac function, represents shifting the spectrum to the base frequency band;

[0017] Variational solution: Introduce the Lagrange multiplication operator λ and the penalty coefficient α, convert the constrained variational problem into an unconstrained problem, and then use the alternating direction method of multipliers (ADMM) to continuously update the center frequencies w k of each mode, the Lagrange multiplier λ, and the mode function u k to obtain the optimal solution of the constrained model;

[0018] Determine the intrinsic mode function: Through continuous iteration of the sub-step variational solution, when the constraint conditions are met, stop the iteration to obtain each mode u k (t), and by performing the inverse Fourier transform on it, obtain the real part data Pθ (t), each P θ (t) includes a real part X(t) and an imaginary part Y(t). Then, the real part X(t) is the set of IMFs of the EEG signal, and the imaginary part Y(t) is the set of IMFs of the EMG signal.

[0019] Multi-scale processing: Use the moving average method to perform multi-scale analysis on the IMFs, extract the local change characteristics of the EEG-EMG signals, and select to use a sliding window to average the IMFs sequence:

[0020]

[0021] Among them, L(P θ (t), τ) represents the moving average at scale τ, P θ (t) represents each EEG-EMG IMFs sequence, T is the length of the IMFs sequence, and the EEG-EMG IMFs after multi-scale processing are defined as X S and Y S ;

[0022] Transfer spectral entropy calculation: For the EEG-EMG IMFs sequence after multi-scale processing, perform phase space reconstruction processing, and then through two-dimensional Fourier transform, convert the reconstructed sequence into a two-dimensional matrix to perform transfer spectral entropy calculation;

[0023] Upper limb rehabilitation grade assessment: Based on the Spearman correlation coefficient, map the fitting relationship between the significant area value of MSTSE of each patient and the Brunnstrom assessment scale, analyze the FCMC characteristic indexes and change trends that are significantly correlated with the Brunnstrom assessment scale grade, and perform quantitative assessment of the disease grade of stroke patients.

[0024] Furthermore, the steps of collecting signals and performing preprocessing include:

[0025] Combined with the Brunnstrom assessment scale, according to the clinical manifestations of different stroke patients, under the guidance of clinicians, an experimental paradigm for forward flexion of the arm was designed. According to the international 10-20 system, EEGs of stroke patients were collected at leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2 respectively, and sEMGs at the anterior deltoid, middle deltoid, biceps brachii, triceps brachii, and extensor digitorum of the forearm were collected using surface electromyography sensors; preprocessing of the EEG was completed, and the specific operations included: rereferencing, baseline correction, filtering, eye artifact removal, downsampling, and EEG data segmentation; preprocessing operations for the sEMG included: baseline correction, filtering, and EMG data segmentation, etc.

[0026] Furthermore, the steps of the transfer spectral entropy calculation include:

[0027] Perform phase reconstruction on the multi-scale processed sequence, and for X S and Y S expand them into two-dimensional matrices P S and Q S ;

[0028] Transform the time-domain matrices P S and Q S to the frequency domain through two-dimensional Fourier transform to obtain matrices W S and V S ;

[0029] By calculating the transfer spectral entropy of the M-dimensional vectors w S (f) and v S (f), reveal the information transfer between different frequency components.

[0030] Furthermore, the steps of the variational solution include:

[0031] Introduce the Lagrange multiplier operator λ and the penalty coefficient α:

[0032]

[0033] In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the Lagrange multiplier to ensure the strictness of the constraint conditions;

[0034] Based on the Plancherel Fourier transform, transform the above formula to the frequency domain for calculation, and continuously update the central frequency w k of each mode, the Lagrange multiplier λ, and the mode function u k , to obtain the optimal solution of the constraint model:

[0035]

[0036] In the formula, ∧ represents the Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

[0037] Combined with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solutions to be protected by the present invention are:

[0038] First, use the binary variational mode decomposition-multi-scale transfer spectral entropy method to characterize the FCMC coupling relationship between the motor cortex of the brain and the corresponding muscle tissues, and combine with the Brunnstrom assessment scale to achieve an objective quantitative assessment of the upper limb rehabilitation level of stroke patients.

[0039] (1) The technical solution of the present invention solves the technical problems that people have always been eager to solve but have never been successful:

[0040] Based on bivariate variational mode decomposition (mBVMD), this paper synchronously decomposes electroencephalomyographic (EEG-EMG) signals to obtain EEG-EMG signals in two identical frequency bands, realizes the co-frequency coupling analysis of EEG-EMG signals, and explores the neuromuscular coupling characteristics. Compared with the VMD algorithm, BVMD extends one-dimensional signals to two dimensions, ensures the matching of intrinsic mode function components in terms of quantity and scale, and solves the modal calibration problem between signals. At the same time, based on the phase space reconstruction technology, the frequency domain information of EEG-EMG signals is extracted, and the transfer spectral entropy of stroke patients at different levels is calculated, thereby realizing the quantitative and objective assessment of the disease levels of stroke patients.

[0041] (2) The technical solution of the present invention overcomes the technical prejudice:

[0042] In the clinical assessment process of the disease level of stroke, clinicians complete the assessment of the disease level of stroke patients based on the rating scale and clinical subjective experience, which has certain subjectivity. The present invention collects the synchronous EEG-EMG signals of stroke patients during the execution of rehabilitation movements, explores the coupling degree between the signals, maps the relationship between different assessment levels and signal coupling, and realizes the assessment of the disease level of stroke patients. To a certain extent, it overcomes the technical prejudice of subjectively assessing the disease level of stroke patients in the current clinic. Description of the Drawings

[0043] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on the structures shown in these drawings.

[0044] Figure 1 It is the flowchart of the upper limb rehabilitation assessment method for stroke patients based on EEG-EMG signals described in the present invention;

[0045] Figure 2 It is a schematic diagram of the experimental paradigm;

[0046] Figure 3 It is the execution flowchart of the experimental paradigm;

[0047] Figure 4 It is the schematic diagram of the EEG leads used in the present invention;

[0048] Figure 5 It is the schematic diagram of the EMG tissue used in the present invention;

[0049] Figure 6 It is the waveform diagram of the IMFs decomposed from the EEG-EMG complex signal;

[0050] Figure 7 It is the heat map of the MSTSE of the cerebral myoelectric signal;

[0051] Figure 8 It is the mapping relationship diagram between BVMD-MSTSE and the Brunnstrom assessment scale for stroke patients. Specific implementation mode

[0052] The present invention proposes an upper limb rehabilitation assessment method for stroke patients based on cerebral myoelectric signals, aiming to design a rehabilitation assessment method for objectively and quantitatively assessing the upper limb rehabilitation level of stroke patients.

[0053] The following will illustrate the upper limb rehabilitation assessment method for stroke patients based on cerebral myoelectric signals proposed by the present invention in specific embodiments:

[0054] In the technical solution of this embodiment, an upper limb rehabilitation assessment method for stroke patients based on cerebral myoelectric signals includes the following steps:

[0055] Step 1: Signal acquisition and preprocessing;

[0056] Step 2: Decomposition of cerebral myoelectric signals:

[0057] Construct a complex signal, combine the electroencephalogram signal x(t) and the electromyogram signal y(t) into a complex form, defined as: z(t) = x(t) + jy(t), where j is the imaginary unit;

[0058] Direction selection, select a group of directions to project z(t), select multiple different directions θ d , where θ d is the uniformly distributed angle on the unit circle. In the case of selecting D directions, each direction θ d can be expressed as: d = 0, 1, 2,..., D - 1;

[0059] Projection calculation, project the complex signal z(t) onto a certain direction θ, and calculate the projection p θ (t) of z(t) in this direction. p θ (t) = Re(z(t)e -jθ ), where Re represents taking the real part, and θ is the angle of the projection direction;

[0060] Parameter initialization, initialize the parameters required for VMD, including the number of modes K, the central frequency w k of each mode, the Lagrange multiplier λ, the mode function u k , the penalty parameter α, and the convergence tolerance tol;

[0061] Define the objective function, requiring each mode u k(t) has a finite bandwidth and a center frequency, with the constraint that the sum of each mode is equal to the input signal p θ (t), and the sum of the estimated bandwidths of each mode is minimized. The constrained variational model is as follows:

[0062]

[0063] s.t. ∑ k u k = p θ (t);

[0064] In the formula, {w k} represents the frequency center of each IMF component, {u k} are the K finite-bandwidth IMF components obtained after decomposition, is to take the partial derivative, * is the convolution operator, δ(t) is the Dirac function, represents shifting the spectrum to the base frequency band;

[0065] For variational solution, introduce the Lagrange multiplier operator λ and the penalty coefficient α, convert the constrained variational problem into an unconstrained problem, and then use the alternating direction method of multipliers (ADMM) to continuously update the center frequency w k of each mode, the Lagrange multiplier λ, and the mode function u k to obtain the optimal solution of the constrained model;

[0066] Determine the intrinsic mode function. Through continuous iteration of the sub-step variational solution, when the constraint conditions are met, stop the iteration to obtain each mode u k (t). By performing the inverse Fourier transform on it, obtain the real part data P θ (t). Each P θ (t) contains a real part X(t) and an imaginary part Y(t). Then the real part X(t) is the IMF set of the electroencephalogram signal, and the imaginary part Y(t) is the IMF set of the electromyogram signal.

[0067] Step 3: Multi-scale processing. Use the moving average method for multi-scale analysis of the IMFs to extract the local change characteristics of the electroencephalogram and electromyogram signals. Select to use the sliding window to average the IMFs sequence:

[0068]

[0069] where, L(P θ (t), τ) represents the moving average at scale τ, P θ (t) represents each electroencephalogram and electromyogram IMFs sequence, T is the length of the IMFs sequence, and define the electroencephalogram and electromyogram IMFs after multi-scale processing as X S and Y S ;

[0070] Step 4: Transfer spectral entropy calculation. For the brain-muscle EMG IMFs sequence after multi-scale processing, perform phase space reconstruction processing, and then through two-dimensional Fourier transform, convert the reconstructed sequence into a two-dimensional matrix to calculate the transfer spectral entropy;

[0071] Step 5: Upper limb rehabilitation grade assessment. Based on the Spearman correlation coefficient, map the fitting relationship between the MSTSE significant area value of each patient and the Brunnstrom assessment scale, analyze the FCMC characteristic indexes and change trends that are significantly correlated with the Brunnstrom assessment scale grade, and quantitatively assess the disease grade of stroke patients.

[0072] Furthermore, the steps of collecting signals and performing preprocessing include:

[0073] Combined with the Brunnstrom assessment scale, according to the clinical manifestations of different stroke patients, an experimental paradigm of arm flexion is designed under the guidance of a clinician. According to the international 10-20 system, EEGs of stroke patients are collected at leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2 respectively, and sEMGs at the anterior deltoid, middle deltoid, biceps brachii, triceps brachii, and extensor digitorum of the forearm are collected using surface electromyography sensors; preprocessing of EEG is completed, and the specific operations include: rereferencing, baseline correction, filtering, eye electrooculogram removal, downsampling, and electroencephalogram data segmentation; preprocessing operations for sEMG include: baseline correction, filtering, electromyogram data segmentation, etc.

[0074] Furthermore, the steps of calculating the transfer spectral entropy include:

[0075] Perform phase reconstruction processing on the multi-scale processed sequence, and expand X S and Y S into two-dimensional matrices P S and Q S ;

[0076] Convert the time-domain matrices P S and Q S to the frequency domain through two-dimensional Fourier transform to obtain matrices W S and V S ;

[0077] By calculating the transfer spectral entropy of the M-dimensional vectors w S (f) and v S (f), reveal the information transfer between different frequency components.

[0078] Furthermore, the steps of variational solution include:

[0079] Introduce the Lagrange multiplier operator λ and the penalty coefficient α:

[0080]

[0081] In the formula, α is a penalty function to ensure the reconstruction accuracy of the signal, and λ is a Lagrange multiplier to ensure the strictness of the constraint condition;

[0082] Based on the Plancherel Fourier transform, the above formula is transformed into the frequency domain for calculation, and the central frequency w of each mode, k the Lagrange multiplier λ, and the mode function u k are continuously updated to obtain the optimal solution of the constraint model:

[0083]

[0084] In the formula, ∧ represents the Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

[0085] Example 1:

[0086] A method for upper limb rehabilitation assessment of stroke patients based on electroencephalogram and electromyogram signals, and the specific process is as Figure 1 shown, including the following steps:

[0087] Step 1: Signal acquisition and preprocessing. Under the guidance of a clinician, setting the arm flexion forward as an experimental paradigm in combination with the Brunnstrom assessment scale, collecting synchronous electroencephalogram signals and surface electromyogram signals of several stroke patients when performing arm flexion forward, and ensuring that there are at least 5 subjects in each level of Brunnstrom among several stroke subjects.

[0088] In the experimental paradigm, there are two actions: arm flexion forward and rest, as Figure 2 shown. When the subject performs arm flexion forward, a prompter is placed in front of the subject, and the subject is given picture and sound prompts to prompt the subject to perform the corresponding action. During the signal acquisition process, the subject needs to sit upright on the chair, with the arm hanging naturally, keeping the body stable, and reducing unnecessary movements, including but not limited to: head turning, swallowing, blinking, etc.; as Figure 3The specific action execution process is shown as follows. In the preparation stage from 0 to 3 s, the subject is given a 1-s voice prompt, and at the same time, a "+" is displayed on the prompter to focus the subject's attention. During this period, the subject's arm hangs naturally; in the action execution stage, the subject needs to slowly flex the arm forward until it is parallel to the ground and maintain stability. The time for the action execution stage is 8 s; in the rest stage, the subject needs to slowly lower the arm until it hangs naturally and is perpendicular to the ground. The time for the rest stage is 12 s. Each subject needs to complete 3 sessions of actions. Each session of actions includes 10 arm flexion actions and 10 rest actions. After the subject completes 1 session of actions, a 10-minute rest time will be given to prevent the subject from experiencing mental and muscle fatigue.

[0089] During the acquisition process of the EEG, in strict accordance with the international 10-20 system, the subject needs to complete scalp cleaning in advance, correctly wear the EEG cap, and select the central area lead Cz as the reference point. To avoid distortion of the collected EEG signals, the impedance between the electrode and the scalp must be less than 5 kΩ, and then the EEG signals are collected, as Figure 4 shown, the EEGs at the Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2 leads are collected respectively.

[0090] During the acquisition process of the sEMG, the skin of the subject's upper arm is cleaned with alcohol, and if necessary, the local skin of the subject is polished with a scrub to reduce the skin impedance. As Figure 5 shown, the sEMGs at the anterior deltoid, middle deltoid, biceps brachii, triceps brachii, and extensor digitorum of the forearm are collected respectively.

[0091] The preprocessing operations of the EEG described above include but are not limited to rereferencing, baseline correction, filtering, electrooculogram removal, downsampling, and EEG data segmentation, etc. In the rereferencing operation, A1 and A2 are selected as the new reference points; in the baseline correction operation, the baseline drift of the EEG is reduced by the demeaning method; in the filtering operation shown, a band-pass filter of 0.1 - 200 Hz and a notch filter of 50 Hz are completed for the EEG; in the electrooculogram removal operation, the independent component analysis (ICA) algorithm is used to remove the electrooculogram artifacts of the subject; in the EEG signal downsampling operation, the sampling rate of the EEG signal needs to be reduced to the same as that of the electromyogram signal; in the EEG data segmentation operation, the time range of the action execution stage is determined based on the marks in the EEG acquisition system, and the corresponding data segments are extracted.

[0092] The preprocessing operations of the sEMG include but are not limited to baseline correction, filtering, and EMG data segmentation. In the shown baseline correction operation, the baseline drift of the sEMG is reduced by the mean method; in the filtering operation, band-pass filtering of 0.1 - 200 Hz and notch processing of 50 Hz are performed on the sEMG; in the EMG data segmentation operation, the time range of the action execution phase is determined based on the marks in the EMG acquisition system, and the corresponding data segments are extracted.

[0093] Step 2: Decomposition of the brain-EMG signal. In this embodiment, the binary variational mode decomposition (BVMD) is used to realize the synchronous decomposition of the brain-EMG signal, which specifically includes the following contents:

[0094] 2.1 Construction of a complex signal: The preprocessed EEG signal x(t) and the surface EMG signal y(t) are combined into a complex signal z(t) = x(t) + jy(t), where j is the imaginary unit;

[0095] 2.2 Direction selection: Select a set of directions to project z(y). Generally, multiple different directions θ d , where θ d is the uniformly distributed angle on the unit circle. In the case of selecting D directions, each direction θ d can be expressed as: d = 0, 1, 2, …, D - 1, and the value range of D is set between 10 and 20;

[0096] 2.3 Projection calculation: Project the complex signal z(t) onto a certain direction θ, and calculate the projection signal p θ (t) of z(t) in this direction. p θ (t) = Re(z(t)e -jθ ), where Re represents taking the real part, and θ is the angle of the projection direction;

[0097] 2.4 Parameter initialization: Initialize the parameters required for VMD, including the number of modes K, the central frequency w k of each mode, the Lagrange multiplier λ, the mode function u k , the penalty parameter α, and the convergence tolerance tol;

[0098] 2.5 Definition of the objective function: It is required that each mode u k (t) has a finite bandwidth and central frequency, and the constraint condition is that the sum of all modes is equal to the input signal p θ (t), and the sum of the estimated bandwidths of each mode is minimized. The constrained variational model is as follows:

[0099]

[0100] s.t. ∑k u k = p θ (t);

[0101] In the formula, {w k} represents the frequency center of each IMF component, and {u k} are the K finite-bandwidth IMF components obtained after decomposition. ∂ / ∂ is for partial derivative, * is the convolution operator, δ(t) is the Dirac function. represents shifting the spectrum to the base frequency band;

[0102] 2.6 Variational solution. To solve the above constrained optimization problem, introduce the Lagrange multiplier λ and the penalty coefficient α to convert the constrained variational problem into an unconstrained problem:

[0103]

[0104] In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the Lagrange multiplier to ensure the strictness of the constraint condition;

[0105] Then, further solve the variational problem through the alternating direction method of multipliers (ADMM). Based on the Plancherel Fourier transform, transform the above formula to the frequency domain for calculation, and continuously update the central frequency w k of each mode, the Lagrange multiplier λ, and the mode function u k to obtain the optimal solution of the constrained model:

[0106]

[0107] In the formula, ∧ represents the Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

[0108] 2.7 Determine the intrinsic mode function. Set the constraint conditions:

[0109]

[0110] Through continuous iteration of sub-step 2.6, when the constraint conditions are met, stop the iteration to obtain each mode u k (t). By performing the inverse Fourier transform on it, obtain the real part data P θ (t). Each P θ (t) contains the real part X(t) and the imaginary part Y(t). Then, the real part X(t) is the IMF set of the EEG signal, and the imaginary part Y(t) is the IMF set of the EMG signal.

[0111] Step 3: Multi-scale processing. Use the moving average method to perform multi-scale analysis on the IMF and extract the local change characteristics of the EEG-EMG signal:

[0112]

[0113] Among them, L(P θ (t), τ) represents the moving average at scale τ, and P θ (t) represents each electroencephalomyogram (EEG) IMFs sequence. T is the length of the IMFs sequence. Define the EEG IMFs after multi-scale processing as X S and Y S .

[0114]

[0115] Step 4: Calculate the transfer spectral entropy. For the EEG IMFs sequence after multi-scale processing, perform phase space reconstruction processing, and then through two-dimensional Fourier transform, convert the reconstructed sequence into a two-dimensional matrix to calculate the transfer spectral entropy.

[0116] 4.1 To better capture the complex dynamic characteristics and frequency domain characteristics in the signal, first perform phase reconstruction processing on the multi-scale processed sequence, and expand X S and Y S into two-dimensional matrices P S and Q S :

[0117]

[0118] In the formula, σ = 1, 2, …, N - (M - 1)Δt, where M is the embedding dimension, usually set to 15, and Δt is the embedding delay, generally set to 2.

[0119] 4.2 Convert the time-domain matrices P S and Q S to the frequency domain through two-dimensional Fourier transform to obtain matrices W S and V S :

[0120]

[0121] In the formula, the frequency resolution is Δf = fs / N, fs is the sampling frequency of the EEG signal, and N = T - τ + 1.

[0122] 4.3 Reveal the information transfer between different frequency components by calculating the transfer spectral entropy of the M-dimensional vectors w S (f) and v S (f):

[0123]

[0124] Step 5: Evaluate the upper limb rehabilitation level.

[0125] In the upper limb rehabilitation level assessment, the consistency relationship between the transfer spectral entropy and the clinical assessment scale is explored based on the Spearman method. The calculation process of the Spearman correlation coefficient is as follows:

[0126]

[0127] where ρ is the Spearman correlation coefficient, and h i is the rank difference of each group of data. In the present invention, the result of each group of data is recorded as the correlation with the Brunnstrom level. num is the total number of observed values. After calculating the MSTSE significant area values between each muscle tissue and the brain respectively, the mean value of ρ is obtained mean to realize the characterization of the fitting relationship between the MSTSE significant area value of each patient and the Brunnstrom assessment scale, analyze the FCMC characteristic indexes and change trends that are significantly correlated with the Brunnstrom assessment scale level, and complete the quantitative assessment of the disease level of stroke patients.

[0128] Step 6: Instance test

[0129] The present invention collects the electroencephalomyogram signals of 6 stroke patients when performing the experimental paradigm. All patients have right upper limb hemiplegia without cognitive impairment and no history of any mental diseases. All patients are informed of the experimental protocol before the experiment and have never had similar experimental experiences. Among all stroke patients, subject 1 corresponds to a stage 6 patient in the Brunnstrom upper limb assessment scale, subjects 2 and 5 correspond to stage 5 patients in the Brunnstrom upper limb assessment scale, subjects 3 and 6 correspond to stage 4 patients in the Brunnstrom upper limb assessment scale, and subject 4 corresponds to a stage 3 patient in the Brunnstrom upper limb assessment scale.

[0130] The preprocessed action segment electroencephalomyogram signals are combined into complex signals, and the decomposition of the complex signals is realized based on BVMD. Since the most significant brain region controlling human limb movement is located in the contralateral motor area of the brain, the present invention selects the electroencephalogram at the C3 lead of the motor area and the electromyogram signals for coupling characteristic analysis. As Figure 6 shown is the waveform diagram of IMFs after decomposing a certain complex signal, and 5 groups of IMFs signal waveform diagrams are obtained.

[0131] After using the moving average method to complete the multi-scale processing of the IMFs signals, the coupling strength between the multi-scale signals is calculated based on the transfer spectral entropy. As Figure 7The figure shows the distribution schematic diagram of calculating the average value of transfer spectral entropy between the brain-muscle electrical signals of 6 subjects. The vertical coordinate represents each EEG IMF, and the horizontal coordinate successively represents the IMFs of the anterior deltoid muscle, middle deltoid muscle, biceps brachii, triceps brachii, and extensor digitorum of the forearm. It can be seen from the figure that the coupling strength of MSTSE of subject 1 is the strongest, and that of subject 4 is the weakest. Since the Brunnstrom stage of subject 1 is stage 6, when performing the paradigm movement, its muscle contraction state is similar to that of normal people. The coupling strength between the anterior deltoid muscle, middle deltoid muscle and the brain is the strongest, the coupling strength between the biceps brachii, triceps brachii and the brain is the second, and the coupling strength between the extensor digitorum of the forearm and the brain is the lowest; the Brunnstrom stage of subject 4 is stage 3, and when performing the paradigm movement, there are phenomena of shoulder adduction and elbow flexion, which is the most different from the muscle contraction state of normal people. The coupling strength between the biceps brachii, triceps brachii and the brain is the strongest, the coupling strength between the anterior deltoid muscle, middle deltoid muscle and the brain is the second, and the coupling strength between the extensor digitorum of the forearm and the brain is the lowest.

[0132] By calculating the significant area values of MSTSE of stroke patients at each level and mapping the correlation between each MSTSE and the Brunnstrom assessment scale based on Spearman correlation analysis, as Figure 8 shown in the mapping relationship diagram. It can be seen from the figure that there is a significant correlation between the significant area value of MSTSE between brain-muscle electrical signals and the Brunnstrom assessment scale, and there is a positive correlation, indicating that the FCMC characteristics calculated based on BVMD-MSTSE can be used to evaluate the rehabilitation status of stroke patients.

[0133] As mentioned above, it is only the preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed by the present invention should be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. A method for evaluating upper limb rehabilitation of stroke patients based on brain electromyographic signals, characterized in that: The following steps are involved: Signal acquisition and preprocessing; Brain and myoelectric signal decomposition: Construct a complex signal, combine the EEG signal x(t) and the EMG signal y(t) into a complex form, defined as: z(t) = x(t) + jy(t), where j is an imaginary unit; Direction selection, select a set of directions to project the complex signal z(t), select multiple different directions θ d , where θ d is a uniformly distributed angle on the unit circle, select D directions, each direction θ d It can be expressed as: Projection calculation, project the complex signal z(t) to a certain direction θ, and calculate the projection p of z(t) in this direction θ (t), p θ (t) = Re (z (t) e -jθ ), where Re represents the real part and θ is the angle of the projection direction; Parameter initialization, initialization of the parameters required for VMD, including the number of modes K, the center frequency w of each mode k , Lagrange multiplier λ, modal function u k , penalty parameter α and convergence tolerance tol; Define the objective function, requiring each mode u k (t) has a finite bandwidth and center frequency, with the constraint that the sum of the modes is equal to the input signal p θ (t), the sum of the estimated bandwidths of each mode is minimized, and the constrained variational model is as follows: In the formula, {w k } represents the frequency center of each intrinsic mode function (IMFs), {u k } are the K finite bandwidth IMFs components obtained after decomposition, To find the partial derivative, * is the convolution operator, δ(t) is the Dirac function, Indicates that the spectrum is shifted to the baseband; Variational solution, introduce Lagrange multiplication operator λ and penalty coefficient α, alternating direction multipliers, and continuously update the center frequency w of each mode k , Lagrange multiplier λ and modal function u k , used to obtain the optimal solution of the constraint model; Determine the intrinsic mode function, set the constraints through continuous iteration of step variational solution: When the constraints are met, stop the iteration and get each mode u k (t); perform inverse Fourier transform to obtain the real part data P θ (t), each P θ (t) contains the real part X(t) and the imaginary part Y(t), where the real part X(t) is the IMFs set of the EEG signal and the imaginary part Y(t) is the IMFs set of the EMG signal; Multi-scale processing: Use the moving average method to analyze IMFs at multiple scales, extract the local change characteristics of brain and myoelectric signals, and choose to use a sliding window to average the IMFs sequence: Among them, L(P θ (t),τ) represents the moving average under scale τ, P θ (t) represents each brain and myoelectric IMFs sequence, T is the length of the IMFs sequence, and the brain and myoelectric IMFs after multi-scale processing are defined as X S With Y S ; Transfer spectral entropy calculation: After completing multi-scale processing, the EMG IMFs sequence is reconstructed in phase space, and then converted into a two-dimensional matrix through two-dimensional Fourier transform to calculate the transfer spectral entropy; Upper limb rehabilitation grade assessment: Based on the Spearman correlation coefficient, the fitting relationship between the MSTSE significant area value of each patient and the Brunnstrom rating scale was mapped, and the FCMC characteristic indicators and change trends that were significantly correlated with the Brunnstrom rating scale grade were analyzed to quantitatively assess the disease grade of stroke patients; The step of calculating the transfer spectrum entropy comprises: Perform phase reconstruction on the multi-scale processed sequence and convert X S With Y S Expanded to a two-dimensional matrix P S With Q S ; The time domain matrix P S With Q S Through the two-dimensional Fourier transform to the frequency domain, we get the matrix W S and V S ; By calculating the M-dimensional vector w S (f) and v S (f) The transfer spectral entropy is used to reveal the information transfer between different frequency components.

2. The method for evaluating upper limb rehabilitation of stroke patients based on brain myoelectric signals according to claim 1, characterized in that: The signal acquisition and preprocessing steps include: Combined with the Brunnstrom rating scale, according to the clinical manifestations of different stroke patients, an arm flexion experimental paradigm was designed under the guidance of clinicians. According to the international 10-20 system, EEG of stroke patients was collected at leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2, and surface electromyography sensors were used to collect sEMG at the anterior bundle of the deltoid muscle, middle bundle of the deltoid muscle, biceps brachii, triceps brachii, and extensor digitorum of the forearm. The EEG was preprocessed, including re-referencing, baseline correction, filtering, removal of oculoculography, downsampling, and EEG data segmentation; the sEMG was preprocessed, including baseline correction, filtering, and EMG data segmentation.

3. The method for evaluating upper limb rehabilitation of stroke patients based on brain myoelectric signals according to claim 1, characterized in that: The variational solution step includes: Introduce the Lagrange multiplication operator λ and penalty coefficient α: In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the lagrange multiplier; Based on Plancherel Fourier transform, the above formula is converted to the frequency domain for calculation, and the center frequency w of each mode is continuously updated. k , Lagrange multiplier λ and modal function u k , and obtain the optimal solution of the constraint model: Where ∧ represents Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

Images