A method for determining a motor function level assessment model

Abstract

The disclosed embodiment provides a method for determining a motor function level assessment model, which relates to the field of biological signal detection technology. The specific implementation of the method includes: collecting oxygenated hemoglobin concentration change signals, acceleration signals and angular acceleration signals, and constructing a motor level label; preprocessing various signals to obtain brain network characteristic indicators and kinematic characteristic indicators; calculating the correlation between characteristic indicators and scale evaluation parameters, and screening target characteristic indicators; using the target characteristic indicator or the oxygenated hemoglobin concentration change signal of the corresponding channel as input, and the motor level label as output for training to obtain a motor function level assessment model. This implementation can fully consider the synergy between brain regions and the smoothness and stability of limb movements, and can comprehensively reflect the user's motor function from the perspective of brain-limb synergy, and can provide a high reference value for rehabilitation evaluation.

Description

Technical Field

[0001] The present disclosure relates to the technical field of biological signal detection, and in particular to a method for determining a motor function level assessment model. Background Art

[0002] With the rapid increase in the number of people with motor dysfunction caused by neurological diseases such as stroke and cerebral palsy, the accuracy of motor function assessment has been increasingly valued in the field of assessment as a necessary means to understand the functional level of patients, formulate rehabilitation training programs, and evaluate rehabilitation effects. At present, clinical motor function assessment mainly relies on various motor function assessment scales, and the assessors score the assessment parameters of each scale according to the execution of the movements of the people with motor dysfunction. This not only requires a high cost of time, manpower and material resources, but also has large subjective differences among assessors, resulting in poor assessment accuracy.

[0003] The movement disorders of people with movement dysfunction are usually caused by brain damage or maldevelopment. Therefore, the level of motor function of this group of people is actually closely related to the recovery of their brain function. However, the existing assessment scales do not evaluate the recovery of brain function, which further leads to low accuracy of motor function assessment. Summary of the invention

[0004] In view of this, the embodiments of the present disclosure provide a motor function level assessment model, which can solve the problem that the existing scale assessment method not only requires high time, manpower and material costs, but also has large subjective differences among assessors and does not consider the actual influencing factors of movement disorders, resulting in poor assessment accuracy.

[0005] To achieve the above objectives, the present disclosure provides a method for determining a motor function level assessment model, comprising:

[0006] Collecting the user's oxygenated hemoglobin concentration change signal, acceleration signal and angular acceleration signal to construct a sports level label of the user;

[0007] Preprocessing the oxygenated hemoglobin concentration change signal, the acceleration signal, and the angular acceleration signal respectively to obtain a brain network characteristic index and a kinematic characteristic index;

[0008] Respectively calculating the correlations between the brain network characteristic index, the kinematic characteristic index and the evaluation parameters of each scale, and selecting target characteristic indexes therefrom;

[0009] The target characteristic index is used as an input of a motor function level assessment model to be trained, or the oxygenated hemoglobin concentration change signal of the channel corresponding to the target characteristic index is used as an input of the motor function level assessment model to be trained, and the sports level label of the user is used as an output of the motor function level assessment model to be trained, and the motor function level assessment model to be trained is iteratively trained;

[0010] According to the training results of the iterative training, the parameter values ​​of the various model parameters are determined to obtain the motor function level assessment model.

[0011] One or more technical solutions provided in the embodiments of the present application, through training a motor function level assessment model, collect the user's oxygenated hemoglobin concentration change signal, acceleration signal and angular acceleration signal, pre-process to obtain target brain network characteristic indicators and target kinematic characteristic indicators, input the motor function level assessment model and then output the user's motor function level, can achieve the technical effect of simultaneously considering the relevant indicators of brain function and motor function, and accurately and comprehensively reflecting the functional level of people with motor dysfunction. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 A flowchart showing a method for determining a motor function level assessment model according to an exemplary embodiment of the present disclosure is shown;

[0013] FIG. 2( a ) shows a schematic diagram of a system for evaluating a motor function level according to an exemplary embodiment of the present disclosure. Figure 1 ;

[0014] FIG2( b ) shows a second schematic diagram of a system for evaluating a motor function level according to an exemplary embodiment of the present disclosure;

[0015] Figure 3 A flowchart showing a method for constructing a sports level label according to an exemplary embodiment of the present disclosure is shown;

[0016] Figure 4 A flowchart showing a method for processing brain network characteristic indicators according to an exemplary embodiment of the present disclosure is shown;

[0017] Figure 5 A flowchart showing a method for determining a functional connectivity matrix between brain regions according to an exemplary embodiment of the present disclosure is shown;

[0018] Figure 6 A flowchart showing a method for calculating brain network topology parameters according to an exemplary embodiment of the present disclosure is shown;

[0019] Figure 7 A flow chart showing a method for determining a dynamic time-varying index according to an exemplary embodiment of the present disclosure is shown;

[0020] Figure 8 A flowchart of a method for determining an inter-brain effect connection matrix according to an exemplary embodiment of the present disclosure is shown;

[0021] Fig. 9 A flowchart showing a method for determining a community internal functional connection matrix according to an exemplary embodiment of the present disclosure is shown;

[0022] Fig.10 A flow chart showing a method for processing kinematic characteristic indicators according to an exemplary embodiment of the present disclosure;

[0023] Fig.11 A flow chart of a method for screening feature indicators according to an exemplary embodiment of the present disclosure is shown. DETAILED DESCRIPTION

[0024] Embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although certain embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be construed as being limited to the embodiments described herein, which are instead provided for a more thorough and complete understanding of the present disclosure. It should be understood that the drawings and embodiments of the present disclosure are only for exemplary purposes and are not intended to limit the scope of protection of the present disclosure.

[0025] It should be understood that the various steps described in the method embodiments of the present disclosure may be performed in different orders and / or in parallel. In addition, the method embodiments may include additional steps and / or omit the steps shown. The scope of the present disclosure is not limited in this respect.

[0026] The term "including" and its variations used in this document are open inclusions, that is, "including but not limited to". The term "based on" means "based at least in part on". The term "in the embodiment of the present disclosure" means "at least one embodiment"; the term "another exemplary embodiment" means "at least one other embodiment". Relevant definitions of other terms will be given in the following description. It should be noted that the concepts of "first", "second", etc. mentioned in the present disclosure are only used to distinguish different devices, modules or units, and are not used to limit the order or interdependence of the functions performed by these devices, modules or units.

[0027] It should be noted that the modifications of "one" and "plurality" mentioned in the present disclosure are illustrative rather than restrictive, and those skilled in the art should understand that unless otherwise clearly indicated in the context, it should be understood as "one or more".

[0028] The names of the messages or information exchanged between multiple devices in the embodiments of the present disclosure are only used for illustrative purposes and are not used to limit the scope of these messages or information.

[0029] HbO2: oxyhemoglobin, also known as oxygenated hemoglobin, is a substance produced by the reversible combination of hemoglobin and oxygen molecules, which can transport oxygen to tissues.

[0030] Jerk entropy: Instantaneous Rate of Change of Entropy is an indicator used to describe the rate of change and degree of uncertainty of a system or process within a given time. The higher the jerk entropy, the greater the rate of change and degree of uncertainty of the system or process; conversely, the smaller the rate of change and degree of uncertainty.

[0031] Isotonic exercises are exercises in which the length of the muscle changes as it contracts, while the tension remains relatively constant.

[0032] Isometric exercise: isometric contraction refers to the movement in which the length of the muscle remains unchanged during contraction, but tension is generated.

[0033] FMA-UE: Fugl-Meyer Assessment Upper Extremity, is a scale for assessing the motor function of the upper limb on the hemiplegic side, and is the upper limb part of the Fugl-Meyer assessment scale.

[0034] Aspects of the present disclosure are described below with reference to the accompanying drawings.

[0035] Figure 1 A flow chart showing a method for determining a motor function level evaluation model according to an exemplary embodiment of the present disclosure is shown. Figure 1 As shown, the method for determining the motor function level assessment model disclosed in the present invention includes the following steps:

[0036] In the embodiment of the present disclosure, the method for determining the motor function level assessment model of the present disclosure is executed by the assessment server of the motor function level assessment system of the present disclosure, such as Figure 2(a) and 2(b)As shown, the motor function level assessment system disclosed in the present invention includes a task guidance module 1, a data acquisition module 2, a data processing and analysis module 3, a motor function assessment module 4 and an assessment result display module 5. Among them, the task guidance module 1, the data processing and analysis module 3, the motor function assessment module 4 and the assessment result display module 5 are integrated in the application software of the main controller 6 for operation, and the guidance paradigm of the task guidance module 1 and the assessment result of the assessment result display module 5 are displayed through the display 7, which is used to demonstrate actions to the user or display assessment results; the data processing and analysis module 3 is used to process various signals to obtain various characteristic indicators, and the motor function assessment module 4 uses the characteristic indicators to assess the user's motor function level. The data acquisition module 2 and the display 7 are connected to the main controller 6 respectively.

[0037] Specifically, the task guidance module 1 can guide the resting state or the moving state task paradigm, and give the subject task instructions in different forms such as sound, text, picture, video, etc., to guide the subject to complete the task required for the assessment. Among them, the moving state task paradigm includes different movement paradigms such as upper limb elbow extension to touch the knee, upper limb elbow flexion to touch the ear, upper limb bilateral rotation movement, lower limb bicycle movement, walking, etc.

[0038] The data acquisition module 2 includes a brain blood oxygen acquisition device 2A based on near-infrared spectroscopy, a motion information acquisition device 2B based on inertial sensing, and a brain blood oxygen and motion information acquisition synchronization trigger device 2C. During use, the brain blood oxygen acquisition device 2A is worn on the human head and is connected to the synchronization trigger device 2C wirelessly. The motion information acquisition device 2B includes an inertial sensor 2B1 and a motion information acquisition device host 2B2. The inertial sensor 2B1 is worn on each limb of the human body and is connected to the motion information acquisition device host 2B2 wirelessly. At the same time, the motion information acquisition device host 2B2 is connected to the synchronization trigger device 2C, and the motion information acquisition device host 2B2 and the synchronization trigger device 2C are also connected to the main controller 6 respectively.

[0039] The data processing and analysis module 3 is used to perform preprocessing such as data format conversion and filtering on the oxygenated hemoglobin concentration change signal, acceleration signal, angular acceleration signal, etc., and calculate brain network characteristic indicators such as the functional connection matrix between brain regions, brain network topological parameters, dynamic time-varying indicators, brain region effect connection matrix and community internal functional connection matrix, as well as kinematic characteristic indicators such as acceleration jerk entropy, isotonic instability, and isometric instability.

[0040] The motor function assessment module 4 includes a motor function grade assessment model, which inputs the brain network characteristic index and kinematic characteristic index calculated by the data processing and analysis module 3 into the motor function grade assessment model, and outputs the motor function grade to achieve motor function assessment.

[0041] The evaluation result display module 5 is used to display an evaluation report based on motor function evaluation, including various brain network characteristic indicators and kinematic characteristic indicators and corresponding motor function levels, etc.

[0042] Step 101 : collecting the user's oxygenated hemoglobin concentration change signal, acceleration signal and angular acceleration signal to construct a motion level label of the user.

[0043] In the disclosed embodiments, for various task paradigms, the resting state refers to the user being in a sitting state; the moving state refers to the user being in different movement paradigms. The movement paradigm can be selectively set according to actual evaluation needs. For example, the movement paradigm includes upper limb elbow extension to touch knee, upper limb elbow flexion to touch ear, upper limb bilateral rotation movement, lower limb cycling movement, walking, etc.

[0044] Furthermore, the cerebral blood oxygen collection device 2A includes multiple near-infrared spectrometer transmitting probes and near-infrared spectrometer receiving probes. A channel is formed between a transmitting probe and a corresponding receiving probe. The cerebral blood oxygen collection device 2A covers multiple channels of the prefrontal lobe and motor area of ​​the brain.

[0045] In the resting state and the task state, the cerebral blood oxygen collection device 2A transmits the oxygenated cerebral hemoglobin concentration change signal (hereinafter referred to as HbO2 signal) of each channel of the user to the main controller 6. According to the Brodmann partitioning of the brain, the bilateral hemispheres can be divided into the prefrontal lobe (i.e. PFC), premotor area (i.e. PMA), supplementary motor area (i.e. SMA), primary motor area (i.e. M1) and sensory motor area (i.e. SMC). Using C to represent the cerebral hemisphere on the same side as the motor limb and I to represent the cerebral hemisphere on the opposite side of the motor limb, the bilateral hemispheres can be divided into 10 brain areas: CPFC, CPMA, CSMA, CM1, CSMC and IPFC, IPMA, ISMA, IM1, ISMC.

[0046] In the motion state, before acquisition, each inertial sensor 2B1 is installed on each body part of the user according to the motion paradigm to be executed; during acquisition, by displaying the motion paradigm to be executed on the display 7, the user is guided to execute the corresponding motion paradigm according to the action of the demonstration character on the display screen, and the inertial sensor 2B1 will collect the acceleration signal, angular acceleration signal, etc. of the user and transmit it to the main controller 6. Each motion paradigm is a series of multi-joint linkage movements, including an action execution phase and an action maintenance phase. The display screen can guide the user to execute each motion paradigm in a cycle of action execution A seconds-action maintenance B seconds-action execution A seconds-action maintenance B seconds. Among them, under each motion paradigm, the number and body parts of the inertial sensor 2B1 can be selectively set as needed. For example, if the motion paradigm is upper limb extension of elbow to touch knee, the inertial sensor 2B1 is installed on the upper arm, forearm, hand, waist and back of the user, and the acceleration signal and angular acceleration signal are collected to characterize the linkage movement of multiple joints such as shoulder joint, elbow joint, wrist joint, etc. under the upper limb extension of elbow to touch knee motion paradigm.

[0047] Furthermore, if Figure 3 As shown, the method for constructing the sports level label of the present disclosure includes the following steps:

[0048] Step 301 : issuing an inspection task, wherein the inspection task includes a plurality of scale assessment parameters of a motor function assessment scale.

[0049] In the disclosed embodiment, the main controller 6 issues inspection tasks according to the various scale evaluation parameters of the motor function evaluation scale, so that the user performs the various inspection tasks.

[0050] Furthermore, the motor function assessment scale is the Fugl-Meyer motor function assessment scale, which includes scale assessment parameters of reflex activity, coordinated movement, wrist, elbow, fingers, etc. of the upper limbs, and scale assessment parameters of reflex activity, coordinated movement, tremor, etc. of the lower limbs. For example, the scale assessment parameters of the reflex activity of the upper limbs include biceps tendon reflex and triceps tendon reflex, and the scale assessment parameters of the flexor coordinated movement of the upper limbs include shoulder elevation, shoulder retraction, abduction, external rotation, elbow flexion, forearm supination, etc.

[0051] Furthermore, the inspection task corresponds to the movement paradigm. When the movement paradigm is upper limb movement, the inspection task is the inspection task of each scale evaluation parameter corresponding to the upper limb; when the movement paradigm is lower limb movement, the inspection task is the inspection task of each scale evaluation parameter corresponding to the lower limb; when the movement paradigm is whole-body movement, the inspection task includes the inspection task of each scale evaluation parameter corresponding to the upper and lower limbs.

[0052] Step 302: When the user performs the inspection task, each scale evaluation parameter of the user is scored.

[0053] In the disclosed embodiment, when the user performs each inspection task, the professional scores each scale evaluation parameter of the user according to the scoring criteria of each scale evaluation parameter.

[0054] Step 303: Determine the user's exercise level label by comprehensively analyzing the scores of the evaluation parameters of the scale.

[0055] In the disclosed embodiment, the user's motor function level is determined by synthesizing the scores of the various scale evaluation parameters of the user as the user's motor level label. The motor function level includes level I, level II and level III. For example, if the movement paradigm is upper limb movement, the level is divided according to the FMA-UE score, with a full score of 66 points. The motor function level of users with a comprehensive score between 0 and 22 points is level I, the motor function level of users with a comprehensive score between 23 and 44 points is level II, and the motor function level of users with a comprehensive score between 45 and 66 points is level III.

[0056] Furthermore, the sports function level is numerically encoded to obtain the user's sports level label, for example, level I is encoded as "00", level II is encoded as "01", and level III is encoded as "11".

[0057] Furthermore, or, users who have a historical motor function assessment scale can be selected, and the historical motor function assessment scale can be directly used to construct the user's sports level label, thereby reducing labor costs while ensuring the accuracy of the assessment.

[0058] In the embodiment of the present disclosure, the method for constructing the motion level label of the present disclosure is used to construct the user's motion level label using a motion function assessment scale, which can be used for model training later to improve the accuracy of motion function assessment.

[0059] Step 102 , preprocessing the oxygenated hemoglobin concentration change signal, the acceleration signal, and the angular acceleration signal, respectively, to obtain brain network characteristic indicators and kinematic characteristic indicators.

[0060] In the disclosed embodiment, after receiving the oxygenated hemoglobin concentration change signal sent by the cerebral blood oxygen collection device 2A and the acceleration signal and angular acceleration signal sent by the inertial sensor 2B1, the main controller 6 pre-processes each signal to obtain the user's brain network characteristic indicators and kinematic characteristic indicators.

[0061] Furthermore, during preprocessing, one of the oxygenated hemoglobin concentration change signals in a resting state or the oxygenated hemoglobin concentration change signals in a moving state is preferentially selected to obtain brain network characteristic indicators.

[0062] Alternatively, during preprocessing, the oxygenated hemoglobin concentration change signal in a resting state and the oxygenated hemoglobin concentration change signal in a moving state may be processed simultaneously to obtain brain network characteristic indicators.

[0063] Furthermore, if Figure 4 As shown, the method for processing brain network characteristic indicators disclosed in the present invention includes the following steps:

[0064] In the disclosed embodiments, the characteristic indicators of the brain network include the functional connectivity matrix between brain regions, the topological parameters of the brain network, the effect connectivity matrix between brain regions, and the functional connectivity matrix within the brain network community, wherein the functional connectivity matrix between brain regions characterizes the functional connectivity relationship between the brain regions of interest, and is determined by calculating the phase-locking coefficient between the HbO2 signals of the two channels; the topological parameters of the brain network include the clustering coefficient, the global efficiency, the local efficiency, the transmission efficiency, the average path length, and the network diameter, which are determined by the functional connectivity values ​​between brain regions; the effect connectivity matrix between brain regions is obtained by quantifying the causal relationship between the channel pairs to obtain the effect connectivity matrix between brain regions composed of the causal values ​​of the brain regions; the functional connectivity matrix within the brain network community characterizes the functional connectivity relationship between the communities composed of channels with strong brain connectivity, and is determined by calculating the average value of the phase-locking coefficient between the HbO2 signals of all channel pairs within the two communities.

[0065] Step 401, calculate the phase-locking coefficient between the oxygenated hemoglobin concentration change signals of every two channels as the functional connectivity value between the channels, summarize the functional connectivity values ​​between channels in different brain regions, and determine the functional connectivity matrix between brain regions.

[0066] In the disclosed embodiment, the HbO2 signal refers to the relationship between the amplitude of the brain hemoglobin concentration of each channel in each brain region and time, that is, the HbO2 signal refers to the time series of the brain hemoglobin concentration of each channel. By performing Hilbert transform on the HbO2 signal, calculating the instantaneous phase of each channel, and then calculating the phase-locking coefficient between channels based on the instantaneous phase, and finally taking the average value of the phase-locking coefficients of each multiple channel of interest as the matrix element, the inter-brain functional connection matrix between each brain region is obtained.

[0067] Furthermore, if Figure 5 As shown, the method for determining the functional connection matrix between brain regions disclosed in the present invention includes the following steps:

[0068] Step 501 , performing Hilbert transform on the oxygenated hemoglobin concentration change signal of each channel to obtain the Hilbert value of the oxygenated hemoglobin concentration change signal of each channel.

[0069] In the embodiment of the present disclosure, the HbO2 signal of each channel is subjected to Hilbert transform to obtain the Hilbert value x of the HbO2 signal of each channel.H (t) is shown in the following formula (1):

[0070]

[0071] In the above formula, τ is the variable in the integral, representing the delay time, which is the offset relative to the current time t when calculating the Hilbert transform; x(τ) is the HbO2 signal at the τth moment; PV is the Cauchy principal value; h is the Hilbert transform operator, which represents the Hilbert transform operation of the HbO2 signal x(τ). u is the integral variable, which is used to represent the integral of time t during the integration process.

[0072] Step 502: Calculate the instantaneous phase of the oxygenated hemoglobin concentration change signal of each channel using the Hilbert value.

[0073] In the embodiment of the present disclosure, the Hilbert value x of each channel is used. H (t), calculate the instantaneous phase of each channel As shown in the following formula (2):

[0074]

[0075] In the above formula, x(t) represents the HbO2 signal of each channel at the current time t, that is, the time series of the brain hemoglobin concentration change signal.

[0076] Step 503: Calculate the phase-locking coefficient value between each channel pair according to the instantaneous phase of the oxygenated hemoglobin concentration change signal of each channel as the functional connection value between the channels.

[0077] In the embodiment of the present disclosure, the phase locking coefficient PLV between the channel pair consisting of channel a and channel b is ab As shown in the following formula (3):

[0078]

[0079] In the above formula, PLV ab The larger the value, the more synchronized the phases between channels A and B are, which also means the stronger the functional connectivity between channels A and B is. k is a time variable, representing the kth sampling point in the time series; is the channel a of the first brain region at time t k The instantaneous phase at is the channel b of the second brain region at time t k j is the instantaneous phase at ; j is the index variable of the time series, which is used to iteratively calculate the phase difference at different time points; N is the total number of sampling points or time points.

[0080] Step 504, calculating the average of the functional connectivity values ​​of all channel pairs between the various brain regions as the functional connectivity value between the brain regions, and constructing the functional connectivity matrix between the sample brain regions.

[0081] In the disclosed embodiment, since each brain region has multiple channels, there are also multiple channel pairs between brain regions. Therefore, the functional connectivity values ​​of each channel pair between the two brain regions are calculated respectively, and the functional connectivity values ​​of all channel pairs between the two brain regions are averaged as the functional connectivity value between the two brain regions. By analogy, the functional connectivity values ​​between all brain regions can be obtained, and after combination, the functional connectivity matrix between sample brain regions can be constructed. For example, taking the calculation of the functional connectivity value between brain region CPFC and brain region CPMA as an example, CPFC includes channels CPFC1 and CPFC2, CPMA includes channels CPMA1 and CPMA2, ​​and the channel pairs between CPFC and CPMA include CPFC1 and CPMA1, CPFC2 and CPMA1, CPFC1 and CPMA2, ​​CPFC2 and CPMA2. The phase-locked coefficient of the channel pair CPFC1 and CPMA1 is PLV CPFC1,CPMA1 , the phase-locked coefficient of the channel pair CPFC2 and CPMA1 is PLV CPFC2,CPMA1 , the phase-locked coefficient of the channel pair CPFC1 and CPMA2 is PLV CPFC1,CPMA2 , the phase-locked coefficient of the channel pair CPFC2 and CPMA2 is PLV CPFC2,CPMA2 Therefore, the phase-locking coefficient between CPFC and CPMA is PLV CPF1C ,CPM1A,PLV CPFC2,CPMA1 、PLV CPFC1,CPMA2 and PLV CPFC2,CPMA2 The average value of is:

[0082] (PLV CPFC1,CPMA1 +PLV CPFC2,CPMA1 +PLV CPFC1,CPMA2 +PLV CPFC2,CPMA2 ) / 4.

[0083] Furthermore, the inter-brain functional connectivity matrix R FCBR As shown in Table 1 below:

[0084] Table 1

[0085] CPFC CPMA CSMA CM1 CSMC IPFC IPMA ISMA IM1 ISMC CPFC * * * * * * * * * * CPMA * * * * * * * * * * CSMA * * * * * * * * * * CM1 * * * * * * * * * * CSMC * * * * * * * * * * IPFC * * * * * * * * * * IPMA * * * * * * * * * * ISMA * * * * * * * * * * IM1 * * * * * * * * * * ISMC * * * * * * * * * *

[0086] Step 402, using the inter-brain functional connectivity matrix, calculate brain network topology parameters, including clustering coefficient, global efficiency, local efficiency, transmission efficiency, average path length, and network diameter.

[0087] In the disclosed embodiment, the percentage threshold sparsification method is first used to screen the target phase-locking coefficients at various sparsities from the functional connectivity matrix between brain regions, and then the multiple target phase-locking coefficients are used to calculate the topological parameters of each brain network, and the change curves of each brain network topological parameter are plotted. The change curves are integrated to obtain the clustering coefficient, global efficiency, local efficiency, transmission efficiency, average path length, and network diameter.

[0088] Furthermore, if Figure 6 As shown, the method for calculating brain network topology parameters disclosed in the present invention includes the following steps:

[0089] Step 601, according to the preset sparsity interval and step size, the target phase-locking coefficients satisfying each sparsity in the functional connection matrix between brain regions are screened, and the brain network topology parameters under each sparsity are calculated.

[0090] In the embodiments of the present disclosure, in order to ensure the accuracy of brain network topological parameters and reduce the impact of phase-locked coefficient screening on topological parameter calculation, an appropriate sparsity interval of the percentage threshold sparsification method and an incremental step size of the sparsity within the sparsity interval are set to obtain multiple sparsities corresponding to the sparsity interval. For example, the sparsity interval is [0.2, 0.8] and the step size is 0.01, and the obtained sparsities include 0.2, 0.21, 0.22, ..., 0.79, and 0.8 to remove the impact of weak connections.

[0091] Furthermore, for each sparsity within the sparsity interval, the target phase-locking coefficients less than or equal to the sparsity in the functional connectivity matrix of the brain interval were screened, and the clustering coefficient, global efficiency, local efficiency, transmission efficiency, average path length, and network diameter of the brain network at different sparsities were calculated. Specifically:

[0092] Step 6011, for each sparsity, select a target phase-locking coefficient less than or equal to the sparsity from the inter-brain region functional connectivity matrix.

[0093] In the embodiment of the present disclosure, for example, for a sparsity of 0.2, PLV is screened from the functional connectivity matrix between brain regions. ab Target phase lock coefficient PLV ≤0.2 ab,tar .

[0094] Step 6012, calculating the clustering coefficient under the sparsity to obtain a plurality of clustering coefficient scatter points.

[0095] In the disclosed embodiment, the clustering coefficient C describes the density of all clusters in the brain network, and its value is between 0 and 1. A channel is a node in the brain functional network. The larger the C, the more clustered the nodes in the brain network are, and the higher the local clustering of the brain network. In which, each channel is regarded as a node, and the matrix elements of the functional connection matrix between brain regions are regarded as the edges of the brain network, and the clustering coefficient C of the brain network is shown in the following formula (4):

[0096]

[0097] In the above formula, C i is the clustering coefficient of node i, C is the clustering coefficient of the brain network; E i represents the number of neighboring nodes directly connected to node i, d i It represents the degree of node i, that is, the number of neighbor nodes of node i, reflecting the connection density of node i, and n is the total number of nodes.

[0098] Step 6013, calculating the global efficiency under the sparsity to obtain a plurality of global efficiency scatter points.

[0099] In the disclosed embodiment, the global efficiency refers to the average value of the reciprocal of the shortest path length of all nodes, which characterizes the efficiency of information transmission between nodes in the brain network. The higher the global efficiency, the faster and more efficient the information transmission efficiency between nodes, and the stronger the information integration and coordination capabilities of the brain network. The global efficiency GE of the brain network is shown in the following formula (5):

[0100]

[0101] In the above formula, D ij represents the shortest path length between node i and node j, and n is the total number of nodes.

[0102] Step 6014, calculating the local efficiency under the sparsity to obtain a plurality of local efficiency scatter points.

[0103] In the embodiment of the present disclosure, local efficiency refers to the average efficiency of the subnetwork composed of all neighboring nodes of each node in the brain network, which is used to measure the efficiency of information transmission in the local brain region. The higher the local efficiency, the faster and more efficient the information transmission in the local brain region. The local efficiency LE of the brain network is shown in the following formula (6):

[0104]

[0105] In the above formula, G i represents the subnetwork composed of neighboring nodes of node i; E(G i ) represents the subnetwork G i The efficiency is calculated in the same way as the global efficiency, the only difference is the number of nodes.

[0106] Step 6015, calculating the transfer efficiency under the sparsity to obtain a plurality of transfer efficiency scatter points.

[0107] In the disclosed embodiment, a triple refers to three nodes connected in a brain network, a closed triple (also called a triangle) refers to a triple whose connection shape is a closed triangle, and the transmission efficiency refers to the ratio of the number of closed triples to the total number of triples in the traversal brain network, which reflects the probability that information in the brain network can be quickly transmitted from one node to other nodes. The transmission efficiency T of the brain network is shown in the following formula (7):

[0108]

[0109] In the above formula, Triangle number represents the number of closed triplets, and Triple number represents the total number of triplets.

[0110] Step 6016, calculating the average path length under the sparsity to obtain a plurality of path length scatter points.

[0111] In the embodiment of the present disclosure, the average path length refers to the average value of the distance between any two nodes in the brain network. The average path length L of the brain network is shown in the following formula (8):

[0112]

[0113] Step 6017, calculate the network diameter under the sparsity to obtain multiple network diameter scatter points.

[0114] In the embodiment of the present disclosure, the network diameter refers to the maximum value of the distance between any two nodes in the brain network. The network diameter D of the brain network is shown in the following formula (9):

[0115] D=max{D ij},i=1:n,j=(i-1):n (9)

[0116] Step 602, drawing a variation curve of each brain network topology parameter within the sparsity interval.

[0117] In the embodiment of the present disclosure, each sparsity i in the sparsity interval is used as the horizontal coordinate, and each parameter scatter point F i For example, with 0.2, 0.21, 0.22, ..., 0.79, 0.8, etc. as the horizontal coordinates and clustering coefficient scatter points / global efficiency scatter points / local efficiency scatter points / transmission efficiency scatter points / path length scatter points / network diameter scatter points as the vertical coordinates, the change curves of clustering coefficient, global efficiency, local efficiency, transmission efficiency, path length, and network diameter within the sparsity range are plotted.

[0118] Step 603, respectively integrate the change curves of each brain network topology parameter, calculate the area under each change curve, and obtain the clustering coefficient, global efficiency, local efficiency, transmission efficiency, average path length, and network diameter.

[0119] In the disclosed embodiment, the curve area enclosed by the change curves of each brain network topology parameter is calculated to obtain the clustering coefficient, global efficiency, local efficiency, transmission efficiency, average path length, and network diameter. The calculation method of the curve area enclosed by the change curves of each brain network topology parameter is shown in the following formula (10):

[0120]

[0121] In the above formula, F i is the parameter scatter point of each brain network topology parameter under sparsity i, and M is the parameter value of each brain network topology parameter obtained by area integration.

[0122] Step 403, sparse processing is performed on the functional connection matrix between brain regions to obtain a sparse dynamic functional connection matrix in each time window, matrix elements in different time windows are extracted from the sparse dynamic functional connection matrix to form an element sequence, and the dynamic time-varying index is determined according to the standard deviation of the element sequence.

[0123] In the disclosed embodiment, in order to analyze the dynamic changes of the brain network and better characterize the recovery level of brain function, the functional connection matrix between brain regions is sparsely processed to obtain a dynamic sparse dynamic functional connection matrix; for each channel pair, the matrix elements in the sparse dynamic functional connection matrix under different time windows are extracted to form a dynamic element sequence; the dynamic time-varying index of each channel pair is determined by using each element sequence, and the dynamic time-varying index of the brain network is obtained by combining all channel pairs. Specifically, Figure 7 As shown, the method for determining the dynamic time-varying index disclosed in the present invention includes the following steps:

[0124] Step 701 : Segment the oxygenated cerebral hemoglobin concentration change signal using a sliding time window, and calculate the functional connection matrix between brain regions in each time window.

[0125] In the embodiment of the present disclosure, the window width of the sliding time window can be selectively set as needed. For example, the window width of the sliding time window is 200 sampling points, and the sliding step is 10 sampling points.

[0126] Step 702: sparsely estimate the functional connectivity matrix between brain regions in each of the time windows.

[0127] In the embodiment of the present disclosure, the least absolute shrinkage and selection operator (LASSO) is used to perform sparse processing on the functional connectivity matrix between brain regions in each time window. First, a sparse estimate of the inverse covariance matrix is ​​obtained by minimizing the penalized negative log-likelihood. As shown in the following formula (11):

[0128]

[0129] In the above formula, tr() represents the trace of the matrix; S[n] is the covariance matrix, which is an unbiased estimate of the covariance relationship between variables in the data; C -1 [n] is the estimated value of the inverse covariance matrix; λ is the penalty factor for controlling network sparsity; ||C -1 [n]||1 means C -1 [n] L1 norm; argmin represents the parameter value that makes the expression reach the minimum value, An attempt is made to find the best sparse inverse covariance matrix through an optimization process.

[0130] Step 703: Use a preset selection operator to update the sparse estimation in each of the time windows to determine the sparse dynamic functional connectivity matrix in each of the time windows.

[0131] In the embodiment of the present disclosure, secondly, according to the preset selection operator, the sparse estimation is set or retained. The values ​​of the corresponding position elements of are updated to obtain the sparse dynamic functional connection matrix W[n], as shown in the following formula (12):

[0132]

[0133] In the above formula, w ij [n] is an element of W[n].

[0134] Step 704, for each channel pair corresponding to the matrix element of the sparse dynamic functional connection matrix, extract an element sequence composed of the matrix element values ​​under each time window, use the element sequence to calculate the dynamic time-varying index of each channel pair, and summarize the dynamic time-varying index of each channel pair to obtain the dynamic time-varying index of the brain network.

[0135] In the disclosed embodiment, the dynamic time-varying index of each channel pair is determined by first calculating the deviation between each matrix element value in the element sequence and the average value of the matrix elements in the element sequence; then finding the average value of the sum of the squares of each deviation; and finally taking the square root of the average value.

[0136] Step 404, quantify the causal relationship between each channel pair, obtain the causal value between each channel pair, summarize the causal values ​​of the channel pairs between different brain regions, and construct the inter-brain region effect connection matrix of the brain network.

[0137] In the embodiment of the present disclosure, in order to determine whether there is mutual influence between different channels, the channel signals are preprocessed and input into the multivariate autoregressive model for prediction, and the influence degree between the channels is obtained according to the prediction results, and the causal relationship between each channel pair is obtained, thereby realizing the effect network analysis of the brain network. Specifically, Figure 8 As shown, the method for determining the brain region effect connection matrix disclosed in the present invention includes the following steps:

[0138] Step 801 : Detrending and degrading the oxygenated hemoglobin concentration change signal of each channel to obtain a degraded signal.

[0139] In the disclosed embodiment, the purpose of detrending is to eliminate the influence of the long-term trend of the time series corresponding to the oxygenated hemoglobin concentration change signal on the data analysis and stabilize the variance; the purpose of degradation is to eliminate signal noise and smooth the random noise in the time series.

[0140] Furthermore, the oxygenated hemoglobin concentration change signal can be detrended using the method shown in the following formula (13), and the detrended signal x detrended (t) are as follows:

[0141] x detrended (t) = x(t) - (α + βt) (13)

[0142] In the above formula, x(t) is the time series corresponding to the oxygenated hemoglobin concentration change signal;

[0143] α and β are the intercept and slope respectively obtained by least squares fitting of each sampling point of x(t).

[0144] Step 802 , perform a covariance stationarity test on the degraded signal to determine whether the test result is stationary. If yes, go to step 803 ; if not, go to step 808 .

[0145] In the disclosed embodiment, the verification method of the covariance stationarity test includes the Augmented-Dickey-Fuller test (ADF test for short) and the Kwiatkowski-PhillipsSchmidt-Shin test (KPS test for short), and the user checks whether the time series meets the stationarity condition, that is, the mean and variance of the time series do not change with time to avoid pseudo causality.

[0146] Furthermore, the ADF test and the KPS test test whether the time series has a unit root by constructing a regression model of the time series and testing whether the regression coefficient is equal to 1. When any of the test results of the ADF test and the KPS test is that the time series has a unit root, it is determined that the degradation signal is non-stationary and does not meet the stationarity condition; when the test results of the ADF test and the KPS test are both that the time series has a unit root, it is determined that the degradation signal is stationary and meets the stationarity condition.

[0147] Step 803: Use the information criterion method to select the lag order as the model order of the multivariate autoregressive model.

[0148] In the disclosed embodiments, the information criterion method is generally used to select the optimal model order from time series models of different orders. Specifically, the information criterion method is used to compare the log-likelihood function values ​​and model complexity of time series models using different model orders, thereby determining the optimal model order. After selecting the optimal model order, the optimal model order can be used to construct a multivariate autoregressive model to calculate the causal value between each channel, thereby determining the influence relationship between different channel pairs and different brain regions.

[0149] Furthermore, the present disclosure adopts the Akaike information criterion and the Bayesian information criterion, and selects the maximum model order based on the goodness of fit of the time series model and the complexity of the model for subsequent construction of a multivariate autoregression model (MAR), wherein the maximum model order is the optimal lag order p.

[0150] Step 804, bringing the lag order into the multivariate autoregressive model, for each channel, respectively inputting the time series corresponding to the degradation signal of the channel and the time series corresponding to the degradation signal of the channel pair consisting of the channel and other channels into the multivariate autoregressive model to obtain the first predicted time series and the second predicted time series of the channel.

[0151] In the embodiment of the present disclosure, MAR can use the time series corresponding to the degraded signal of channel a to predict the time series of channel a, and can also use the time series corresponding to the degraded signals of channel a and channel b to predict the time series of channel a at the same time. MAR is shown in the following equations (14) and (15):

[0152]

[0153]

[0154] In the above formula, x(t) a is the time series corresponding to the degradation signal of channel a; x(t) bis the time series corresponding to the degraded signal of channel b; z(t) is the conditional variable, the time series corresponding to the degraded signals of other channels except channel a; α0 is the model intercept term; lag represents the model order; ρ order ,θ order They are x(t) a 、x(t) b , the lag term coefficient of z(t); ε 1t , ε 2t is the error term, usually assumed to be a zero-mean and independent and identically distributed random variable; X(t) a 、X(t) a,b They are the first prediction time series and the second prediction time series respectively.

[0155] Step 805: Perform an F test on the first predicted time series and the second predicted time series according to a preset significance threshold to determine whether there is a causal relationship between the channel pairs.

[0156] In the disclosed embodiment, the Bonferroni correction is a multiple comparison correction method used to control the overall error rate to reduce false positive results. Since it is necessary to calculate the causal values ​​between many channels, it may cause false positive results, so a Bonferroni correction is required; wherein, the significance threshold of the Bonferroni correction can be selectively set as needed, for example, the significance threshold of the Bonferroni correction is 0.05, 0.01, etc., selected according to the strictness of the correction. Under the significance threshold of the Bonferroni correction of 0.01, the residual sum of squares between the first predicted time series and the time series of the channel, and the residual sum of squares between the second predicted time series and the time series of the channel are calculated, and the F check value between the residual sum of squares is calculated to determine whether there is a causal relationship between channel a and channel b.

[0157] Furthermore, the calculation of the F check value is shown in the following formula (16):

[0158]

[0159] In the above formula, E a (t) 2 is X(t) a and x(t) a The residual sum of squares is obtained by a and X(t) a After squaring the difference, summing up all sampling points, we get E a,b (t) 2 is X(t) a,b and x(t) aThe residual sum of squares is ; T is the number of sampling points. Through (T-2p-1) processing, the residual sum of squares can be adjusted based on the total degrees of freedom.

[0160] Furthermore, when the F-check value is less than the significance threshold of 0.01, it is determined that there is a causal relationship between the channel pairs.

[0161] Step 806: quantify the causal relationship between the channel pairs to obtain the causal value between the channel pairs.

[0162] In the embodiment of the present disclosure, the quantified causal value of the causal relationship between channel a and channel b is shown in the following formula (17):

[0163]

[0164] In the above formula, EE a is ε 1t The error variance is

[0165] EE a,b is ε 2t The error variance is

[0166] Furthermore, the MAR model validity test can be performed to check the independence of the model residuals and evaluate the reliability of the model. The model validity test includes the Durbin-Watson test, consistency test and SSR (regression sum of squares) test. The Durbin-Watson test is used to check whether the residuals (error terms) of the model are uncorrelated to avoid excessive deviation between the predicted value and the actual value; the consistency test is used to check whether the model parameters and predictability can remain stable over time; the SSR test is used to evaluate whether the model is overfitting or underfitting.

[0167] Step 807, combining the causal values ​​between the channel pairs between the various brain regions to obtain the brain region effect connection matrix of the brain network.

[0168] In the disclosed embodiment, the causal values ​​between the channel pairs between each brain region are combined to obtain the brain region effect connection matrix of the brain network.

[0169] For any two brain regions of a single subject, by calculating the causal values ​​between the channels of the two brain regions, it is determined whether the proportion of non-zero causal values ​​to all causal values ​​is above 75%. If the proportion of non-zero causal values ​​between the channels of the two brain regions exceeds 75%, it is determined that there is a causal relationship between the two brain regions. For example, LPFC includes channels 1, 2, 3, 4, 5, and 6, and LMC includes channels 10, 11, 12, 13, and 14. The causal values ​​between channels 1, 2, 3, 4, 5, and 6 and channels 10, 11, 12, 13, and 14 are calculated respectively, and 30 causal values ​​are obtained. If more than 23 of the 30 causal values ​​are non-zero, it is considered that there is a causal relationship between LPFC and LMC.

[0170] Furthermore, for subjects at the group level (for example, patients and healthy people are divided, so all patients are considered as one group level, and all healthy people are considered as one group level), when constructing the inter-brain effect connection matrix, in order to further quantify the causal relationship between brain regions, for those with a causal relationship between two brain regions, the average of all causal values ​​between the channels of the two brain regions is taken as the average causal value between the two brain regions; for those without a causal relationship between the two brain regions, the average causal value is 0. In this way, the inter-brain effect connection matrix at the group level can be constructed.

[0171] Step 808 , performing a differential operation on the degraded signal, and then going to step 802 .

[0172] In the embodiment of the present disclosure, when the degraded signal is non-stationary, a differential operation is performed on the time series of the degraded signal, and the degraded signal after the differential operation is subjected to the ADF test and the KPS test again, and so on, until the degraded signal is a stable time series. The number of differential operations is preferably 2, and usually after 2 differential operations, the time series can meet the stationarity condition.

[0173] Step 405 , dividing each of the channels into communities, and using the oxygenated hemoglobin concentration change signal of each channel to calculate the internal functional connection matrix of the community of the brain network.

[0174] In the disclosed embodiment, a channel is used as a node in the brain function network. Each node is first divided into communities, the modularity increment of the node movement is calculated, and then the community is updated based on the modularity increment until the optimal division is obtained.

[0175] Furthermore, if Fig. 9 As shown, the method for determining the internal functional connection matrix of a community disclosed in the present invention includes the following steps:

[0176] Step 901: divide the brain network into multiple initial communities.

[0177] In the disclosed embodiment, each node in the network is regarded as a separate community; for each node, it is assigned to the separate community where each adjacent node is located, that is, each node is combined with the separate community where the adjacent node is located to construct a new community. Thus, the separate community corresponding to each node and the community composed of the node and the adjacent node are respectively used as the initial community to form the initial community structure of the brain network.

[0178] Step 902, calculating the modularity increment of each initial community.

[0179] In the embodiment of the present disclosure, each initial community is taken as the current community, and the modularity increment ΔQ of the current community caused by community movement is shown in the following formula (18):

[0180]

[0181] In the above formula, m represents the sum of the weights of all links in the brain network;

[0182] w i represents the sum of the weights of the links associated with node i in the current community;

[0183] w i,in represents the sum of the weights of the links connecting node i in the current community to other nodes in the current community;

[0184] ∑in represents the sum of the weights of the links in the current community;

[0185] ∑tot represents the sum of the weights on the links connected to the nodes in the current community, that is, it represents the sum of the weights on the links connecting the nodes in the current community with all the nodes inside and outside the current community, including the link weights between the nodes in the current community and the sum of the weights on the links between each node in the current community and the nodes outside the current community.

[0186] Here, a link refers to the connection between two nodes, that is, the connection between every two nodes is counted as a link, and the weight of each link is the phase-locking coefficient between the corresponding node pair.

[0187] Step 903, regard each of the initial communities as a new super node, the edges between the initial communities as the edges between the new super nodes, regard each of the new super nodes, and the new super nodes and adjacent nodes as new initial communities, and calculate the modularity increment of each of the new initial communities.

[0188] In the disclosed embodiment, after calculating the modularity increment, it is necessary to aggregate the nodes in the same community to form a "new super node". This new super node represents the entire community and retains the connection relationship within the community. That is, compression refers to multi-node aggregation, specifically:

[0189] Multiple nodes in the same initial community are compressed and transformed into a new supernode. The edges between nodes in the original community are transformed into edges within the new supernode, and the edges between original communities are transformed into edges between new supernodes. Thus, new initial communities corresponding to each new supernode and new initial communities composed of new supernodes and adjacent nodes are obtained, and the modularity increment of each new initial community is calculated.

[0190] Step 904, determine whether the modularity increment of the new initial community is equal to the previous modularity increment, if yes, go to step 905; if no, go to step 903.

[0191] In the disclosed embodiment, the modularity increment of the new initial community is compared with the modularity increment of the previous initial community to determine whether the modularity increment of the new initial community is equal to the previous modularity increment.

[0192] Step 905: select the community with the largest modularity increment to merge and obtain the optimal community division.

[0193] In the disclosed embodiment, the above steps are repeated until the modularity of each community in the brain network no longer changes, which means that the optimal community division has been achieved. The two communities with the largest increase in modularity are merged into one community, and its maximum value corresponds to the optimal community division, which can be used as the label of the merged community.

[0194] Step 906, calculate the average functional connection values ​​of all channel pairs within each two communities under the optimal community division as the functional connection value within the two communities, combine the functional connection values ​​between all communities, and construct the internal functional connection matrix of the community.

[0195] Furthermore, if Fig.10 The method for processing kinematic characteristic indicators disclosed in the present invention comprises the following steps:

[0196] In the disclosed embodiment, the kinematic characteristic indicators include acceleration jerk entropy, isotonic instability and isometric instability.

[0197] Step 1001, calculate the rate of change of the acceleration signal, perform Fourier transform on the rate of change and calculate the entropy value to obtain the acceleration jerk entropy.

[0198] In the embodiment of the present disclosure, the acceleration jerk entropy is shown in the following equations (19) and (20):

[0199]

[0200] In the above formula, y i is the jerk, which represents the rate of change of the acceleration signal and is the first-order derivative of the acceleration signal with respect to time;

[0201] A i is the jerk y i Fourier transform of ; where p represents the window length of the sliding window, k is the frequency, n is the data length; H represents the acceleration jerk entropy.

[0202] Step 1002: time-integrate the angular acceleration signal during the motion execution phase to obtain the isotonic instability.

[0203] In the embodiment of the present disclosure, the isotonic instability is the angular dispersion of the motion trajectory during the action execution phase, which reflects the degree of non-smoothness of the subject's isotonic motion and the ability to control the motion. The isotonic instability M is integrated with the angular acceleration signal, and the final result reflects the dispersion of the angular acceleration signal, as shown in the following formula (21):

[0204]

[0205] In the above formula, ω z Refers to the amplitude of the angular acceleration signal during the motion execution phase; t1 and t2 are the start and end times of the time series of the angular acceleration signal, respectively.

[0206] Step 1003: time-integrate the angular acceleration signal in the motion maintenance phase to obtain the isometric instability.

[0207] In the embodiment of the present disclosure, the isometric instability is the angular dispersion in the action maintenance phase, which reflects the ability of the subject to maintain limb stability. The isometric instability M is shown in the following formula (22):

[0208]

[0209] In the above formula, ω w Refers to the amplitude of the angular acceleration signal during the motion maintenance phase; t1 and t2 are the start and end times of the time series of the angular acceleration signal, respectively.

[0210] Step 103, respectively calculating the correlations between the brain network characteristic index, the kinematic characteristic index and each scale evaluation parameter, and selecting target characteristic indexes therefrom.

[0211] In the disclosed embodiments, the Pearson correlation coefficients of the functional connectivity matrix between brain regions, clustering coefficient, global efficiency, local efficiency, transfer efficiency, average path length, network diameter, dynamic time-varying index, effect connectivity matrix between brain regions, functional connectivity matrix within communities, acceleration jerk entropy, isotonic instability, isometric instability and the evaluation parameters of each scale are calculated respectively, and the target characteristic indicators that meet the coefficient threshold are screened out.

[0212] Furthermore, if Fig.11 As shown, the screening method of characteristic indicators disclosed in the present invention comprises the following steps:

[0213] Step 1101, respectively calculate the Pearson correlation coefficient between each characteristic index and each scale evaluation parameter.

[0214] Step 1102, statistically analyzing each Pearson correlation coefficient, and screening out target correlation coefficients that are less than a preset significance threshold.

[0215] In the disclosed embodiment, each Pearson correlation coefficient is statistically analyzed, and a hypothesis test (such as a t-test or an F-test) is used to determine statistical significance, and the target correlation coefficient is screened according to a preset significance threshold. For example, if the significance threshold is p<0.05, the Pearson correlation coefficient that is significant at a 95% confidence level is screened as the target correlation coefficient.

[0216] Step 1103: determine the characteristic index corresponding to the target correlation coefficient as the target characteristic index.

[0217] Furthermore, the target characteristic index includes the target brain network characteristic index and / or the target kinematic characteristic index. For example, when performing upper limb extension to touch the knee, upper limb flexion to touch the ear and upper limb bilateral rotation, the target characteristic index includes: ① the brain area functional connection matrix between CSMA-CM1; ② the intra-community functional connection matrix between ISMA, IPMA and IM1; ③ isometric instability.

[0218] Step 104: using the target characteristic index as the input of the motor function level assessment model to be trained, or using the oxygenated hemoglobin concentration change signal corresponding to the target characteristic index as the input of the motor function level assessment model to be trained, and using the user's motor level as the output of the motor function level assessment model to be trained, and iteratively training the motor function level assessment model to be trained.

[0219] In the embodiment of the present disclosure, the algorithm of the motor function level assessment model to be trained may be an LSSVM algorithm or a WTSN-LSTM algorithm, wherein:

[0220] Furthermore, when the algorithm of the motor function level assessment model to be trained is the LSSVM algorithm, the radial basis function is selected as the kernel function, and the radial basis kernel function K(s,s i ) is shown in the following formula (23):

[0221]

[0222] In the above formula, σ is the bandwidth of the Gaussian kernel, σ>0; s i is the feature vector representing the i-th one; s is the feature vector of the new input model.

[0223] The objective function of the model to be trained for motor function level assessment is shown in formula (24):

[0224]

[0225] In the above formula, Y is the model output, γ>0 represents the penalty factor / regularization parameter, baise represents the bias term, α represents the Lagrange multiplier, I is the identity matrix, and Ω is the kernel function K(s,s i ) is a symmetric matrix of order N.

[0226] Furthermore, in the iterative training process of the motor function level assessment model to be trained, the Seagull hyperparameter optimization algorithm is used as a solution algorithm to find the optimal values ​​of the penalty factor γ and the kernel function parameter σ. Specifically:

[0227] First, determine the new position without collision with other seagulls, as shown in equation (25):

[0228] C s (t) = b s ×P s (t) (25)

[0229] In the above formula, P s (t) represents the current position of the seagull corresponding to the tth iteration, C s (t) represents the new conflict-free position of the seagull corresponding to the tth iteration, b s Represents the movement behavior of seagulls in a specified space.

[0230] Secondly, calculate the optimal direction M of the seagull s (t), as shown in the following formula (26):

[0231] M s (t) = B × (P bs (t)-P s (t)),

[0232] In the above formula, P bs(t) represents the optimal position of the seagull, B is a random number; rd represents a random value of [-1,1], which provides an optimization search direction for the seagull algorithm.

[0233] Again, according to the new position and the optimal direction, the final position D of the seagull is determined. s (t), as shown in the following formula (27):

[0234] D s (t)=|C s (t)+M s (t)| (27)

[0235] In the seagull-based optimization algorithm, the attack behavior of seagulls usually refers to the behavior of seagulls competing for food. This attack behavior can be regarded as a high-intensity local search. By simulating this attack behavior, the algorithm can enhance the search intensity near the current optimal solution and improve the ability to find the local optimal solution. Therefore, the local search ability of the algorithm is increased based on the attack behavior of the seagull, and the final position is converted into the specific coordinates of the seagull in the three-dimensional space, as shown in the following formula (28):

[0236]

[0237] In the above formula, r = ue kv .

[0238] Finally, the specific coordinates are brought into the final attack position of the seagull, as shown in the following formula (29):

[0239] P s (t) = D s (t)×x×y×z+P bs (t) (29)

[0240] In the disclosed embodiment, when the model to be trained for motor function level assessment adopts the LSSVM algorithm, the penalty factor γ and the kernel function parameter σ are first randomly assigned values, the output Y of the objective function is calculated, the output Y and the motor function level label are compared, and the loss difference between the two is calculated, so as to reversely modify the penalty factor γ and the kernel function parameter σ, and continuously iterate the training until the optimal penalty factor γ and the kernel function parameter σ are determined.

[0241] Furthermore, when the algorithm of the motor function level assessment model to be trained is the WTSN-LSTM algorithm, it includes:

[0242] The objective function of the model to be trained for motor function level assessment is shown in the following formula (30):

[0243] hp min = arg min r∈R f(r) (30)

[0244] In the above formula, hp min represents the final optimized hyperparameter, and f(r) is the objective function. In the iterative training of the WTSN-LSTM model, f(r) plays the role of evaluating the model performance under different hyperparameter settings, and the optimization of hyperparameters is achieved by minimizing f(r).

[0245] Furthermore, in the iterative training process of the motor function level assessment model to be trained, the Bayesian algorithm is used as a solution algorithm to find the optimal values ​​of hyperparameters such as the number of hidden units, the maximum number of iterations, the size of small batch data, and the initial learning rate in the WTSN-LSTM network. Specifically:

[0246] Assume that the hyperparameter to be optimized is R = {r1, r2, ..., r n}, during the Bayesian optimization iteration, the objective function value of each hyperparameter is calculated separately to obtain the data set D n ={(r1,f(r1)),(r2,f(r2)),……,(r n ,f(r n ))};

[0247] Assume that each data point in the data set obeys a Gaussian distribution, as shown in the following equation (31):

[0248] f(r 1:n )~GP(μ(r 1:n ),∑(r 1:n ,r 1:n )), μ(r)=E(f(r)) (31)

[0249] In the above formula, r 1:n It refers to the set of input points corresponding to sampling point 1 to sampling point n; E() is the expected value of the predicted value at a given point r; GP is a Gaussian distribution.

[0250] Bayes' theorem is shown in the following formula (32):

[0251] P(f(r i+1 )|f(r 1:n ))∝P(f(r 1:n )|f(r i+1 )|)P(f(r i+1 )) (32)

[0252] In the above formula, P(f(r i+1 )) is f(r i+1 ), i=1,2,…,n-1.

[0253] The Excepted Improvement (EI) function is selected as the sampling function, and the next sampling point is found in the data set of the sampled hyperparameters and their corresponding objective function values. The data set of the sampled hyperparameters and their corresponding objective function values ​​is used as historical information, which is the basis for the next sampling and estimation of Bayesian optimization. As the number of iterations increases, the data set will continue to accumulate, thereby improving the model's understanding of the hyperparameter space and helping to find a better hyperparameter combination. The EI acquisition function after adding the trade-off scalar ξ is shown in the following formula (33):

[0254]

[0255] In the above formula, r + is the known hyperparameter value that minimizes the objective function f(r); Φ(Z) is the standard normal cumulative distribution function, which is used to calculate whether the current candidate point reaches r or is better than the current best r + ; σ(r) is the predicted standard deviation of the objective function f(r) at point r; the trade-off scalar ξ is a parameter used in the Bayesian algorithm to balance search and utilization. The selection of sampling points is determined by comparing the expected value of EI(x). If the trade-off scalar ξ is large, it is more inclined to utilize known information; if the trade-off scalar ξ is small, it is more inclined to explore unknown areas. During the iteration process, the trade-off scalar ξ will be dynamically adjusted according to the accumulated historical information to find the global optimal solution.

[0256] In the embodiment of the present disclosure, when the motor function level assessment model to be trained adopts the WTSN-LSTM algorithm, the number of hidden layer neurons, the maximum number of iterations, the batch size, the initial learning rate, the hidden layer weight matrix, the output layer weight matrix, the bias vector, etc. are randomly assigned, and the output f(r) of the objective function is calculated. The output f(r) and the motor function level label are compared, and the loss difference between the two is calculated, so as to reversely modify the hidden layer weight matrix, the output layer weight matrix, and the bias vector, and continuously iterate the training until it is determined that hp min The corresponding hidden layer weight matrix, output layer weight matrix, and bias vector are obtained. Among them, the loss function is loss = ||f(r)-label|| 2 .

[0257] Step 105, determining the parameter value of each model parameter according to the training result of the iterative training, and obtaining the motor function level assessment model.

[0258] In the embodiment of the present disclosure, the motor function level assessment model can be obtained by iteratively training the model to be trained for motor function level assessment and determining the parameter values ​​of the model parameters using different algorithms. During subsequent use, the assessment server 10 can process the subject's HbO2 concentration signal, acceleration signal and angular acceleration signal, screen the target feature indicators, and input them into the motor function level assessment model to determine the subject's motor function level.

[0259] Experimental verification shows that the motor function evaluation model established using the method of this embodiment has greatly improved evaluation accuracy and precision. Compared with using the HbO2 concentration change signal of the whole brain channel as the feature input and adopting the LSTM model, the HbO2 concentration change signal of a specific channel is selected as the feature input in this embodiment. The evaluation accuracy of the motor function evaluation model is increased by 7.5% and the evaluation precision is increased by 8.98%.

[0260] In the embodiments of the present disclosure, the method for determining the motor function level assessment model of the present disclosure fully considers the coordination between brain regions and the smoothness and stability of limb movements, and can comprehensively reflect the user's motor function from the perspective of brain-limb coordination, and can provide a higher reference value for rehabilitation assessment.

Claims

1. A method for determining a motor function level assessment model, characterized in that: include: Collecting the user's oxygenated hemoglobin concentration change signal, acceleration signal and angular acceleration signal to construct a sports level label of the user; Preprocessing the oxygenated hemoglobin concentration change signal, the acceleration signal and the angular acceleration signal respectively to obtain brain network characteristic indicators and kinematic characteristic indicators; the brain network characteristic indicators include the functional connection matrix between brain regions, brain network topological parameters, dynamic time-varying indicators, the effect connection matrix between brain regions, and the internal functional connection matrix of the community; the kinematic characteristic indicators include acceleration jerk entropy, isotonic instability and isometric instability; wherein the oxygenated hemoglobin concentration change signal includes the oxygenated brain hemoglobin concentration change signal of each channel of the user, and the angular acceleration signal includes the angular acceleration signal of the movement execution stage and the angular acceleration signal of the movement maintenance stage; Respectively calculating the correlations between the brain network characteristic index, the kinematic characteristic index and the evaluation parameters of each scale, and selecting target characteristic indexes therefrom; The target characteristic index is used as an input of a motor function level assessment model to be trained, or the oxygenated hemoglobin concentration change signal of the channel corresponding to the target characteristic index is used as an input of the motor function level assessment model to be trained, and the sports level label of the user is used as an output of the motor function level assessment model to be trained, and the motor function level assessment model to be trained is iteratively trained; Determining the parameter values ​​of various model parameters according to the training results of the iterative training to obtain the motor function level assessment model; Wherein, the processing method of the brain network characteristic index includes: Calculating the phase-locking coefficient between the oxygenated hemoglobin concentration change signals of every two channels as the functional connectivity value between the channels, summarizing the functional connectivity values ​​between the channels of different brain regions, and determining the functional connectivity matrix between the brain regions; Calculating the brain network topology parameters using the inter-brain functional connectivity matrix, wherein the brain network topology parameters include clustering coefficient, global efficiency, local efficiency, transmission efficiency, average path length, and network diameter; Performing sparse processing on the functional connection matrix between brain regions to obtain a sparse dynamic functional connection matrix in each time window, extracting matrix elements in different time windows from the sparse dynamic functional connection matrix to form an element sequence, and determining the dynamic time-varying index according to the standard deviation of the element sequence; Quantifying the causal relationship between each channel pair, obtaining the causal value between each channel pair, summarizing the causal values ​​of the channel pairs between different brain regions, and constructing the inter-brain region effect connection matrix; Dividing each of the channels into communities, and calculating the internal functional connection matrix of the community using the oxygenated hemoglobin concentration change signal of each of the channels; The processing method of the kinematic characteristic index includes: Calculating the rate of change of the acceleration signal, performing Fourier transform on the rate of change and calculating an entropy value to obtain the acceleration jerk entropy; Performing time integration on the angular acceleration signal during the motion execution phase to obtain the isotonic instability; The angular acceleration signal in the motion maintenance phase is time-integrated to obtain the isometric instability.

2. The determination method according to claim 1, characterized in that: The motor function level assessment model adopts the LSSVM algorithm.

3. The determination method according to claim 2, characterized in that: The motor function level assessment model adopts the WTSN-LSTM algorithm.

Images