A source imaging inverse problem analysis method and device and electronic equipment

Abstract

The application discloses a kind of source imaging inverse problem analysis method, device and electronic equipment, wherein the method comprises: constructing real brain model and source model;Solving to obtain conduction matrix;According to spatial distance, the dipole in source model is regionalized, and inter-group sparse constraint model is constructed;Dipole intensity generalized total variation sparse constraint model is constructed;Inter-group sparse constraint model and dipole intensity generalized total variation sparse constraint model are fused, and inverse problem solving model is constructed;Based on inverse problem solving model, the activation range and intensity of dipole in source model are determined.The application first regionalizes dipole and considers inter-group sparsity, and describes the structural sparse prior of cortical electrical activity.In addition, the abruptness at the boundary of activated region on cortex and the smoothness within activated region are described by dipole intensity generalized total variation.Finally, the accurate estimation of dipole activation range and intensity is realized by re-weighted solving mode.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of physiological signal processing, and in particular to a source imaging inverse problem analysis method and device and electronic equipment. BACKGROUND

[0002] Electroencephalogram / magnetoencephalogram (EEG / MEG) and other neural electrophysiological signals have the advantages of high time resolution, non-invasiveness, and easy acquisition, making it possible to explore the intracranial brain neural electrical activity on a millisecond time scale through extracranial detection of electromagnetic signals. Therefore, they are a very important brain signal detection / imaging means in clinical and life. However, the spatial resolution of EEG / MEG in detecting neural activity electromagnetic signals is relatively low. In order to accurately obtain the electrical activity on the cerebral cortex, the method of source imaging (M / ESI) is usually used, which includes forward problem modeling and inverse problem solving, that is, first, the brain is discretized and the cortical space neuron electrical activity is modeled as a current dipole (referred to as dipole) model (this model is also called a source model), and then a forward conduction model of the cortical space dipole to the scalp space sensor (i.e., the EEG electrode of the scalp EEG device / the sensor of the brain magnetic) is constructed (the above process is referred to as forward problem modeling), and finally, the cortical space electrical activity is estimated (i.e., inverse problem solving) by using the known M / EEG signal and the constructed forward conduction model. Source imaging combines the subject's head-brain model to improve the spatial resolution of M / EEG to a certain extent, and has important research value in brain-computer interface, epileptic focus localization, etc.

[0003] The existing inverse problem solving cannot accurately estimate the activation range and activation strength of the dipole at the same time, resulting in the inability to estimate the propagation path of the cortical dipole electrical activity in the activation range, and this propagation path can help us find more accurate epileptic focus / decode brain activity. Therefore, it is of great significance to accurately estimate the activation range and activation strength of the dipole and accurately restore the cortical electrical activity for further improving the accuracy of epileptic focus localization / brain activity decoding. SUMMARY

[0004] Therefore, the embodiments of the present application provide a source imaging inverse problem analysis method, device and electronic equipment to solve the defect that the existing inverse problem solving method cannot simultaneously and automatically accurately estimate the activation range and strength of the dipole.

[0005] According to a first aspect, the embodiments of the present application provide a source imaging inverse problem analysis method, comprising:

[0006] constructing a real head-brain model and a source model;

[0007] solving a conduction matrix based on the structure of the electrode / sensor, the real head-brain model, the source model, and the positional relationship between the electrode / sensor and the real head-brain model.

[0008] regionalize the dipoles in the source model according to spatial distance, and construct an inter-group sparse constraint model;

[0009] construct a dipole intensity generalized total variation sparse constraint model;

[0010] fuse the inter-group sparse constraint model and the dipole intensity generalized total variation sparse constraint model, and construct an inverse problem solving model based on the conduction matrix;

[0011] determine the activation range and intensity of the dipoles in the source model based on the inverse problem solving model.

[0012] In some optional embodiments, the inter-group sparse constraint model is wherein, I is the total number of groups of the dipoles in the source model, is the number of dipoles in the i-th group, i is a subset of a set of activation intensities of all dipoles in the i-th group, and is a subset of a set of all elements in a vector of activation intensities of all dipoles in the source model, , , is the activation intensity of the j-th dipole in the i-th group, N s is the activation intensity of the j-th dipole in the i-th group, N s is the number of dipoles in the source model, is a vector of index numbers of the dipoles in the i-th group in all dipoles, i is the number of dipoles in the i-th group. i

[0013] In some optional embodiments, the regionalization of the dipoles in the source model according to spatial distance comprises:

[0014] S201: obtaining the dipoles that have not been regionalized, and selecting one of the dipoles as a seed point;

[0015] S202: selecting dipoles in an n-order neighborhood of the seed point to form a group together with the seed point; wherein, , , is the maximum order, i.e. ​​The dipoles in the order neighborhood of the seed point are all the dipoles in the source model; the dipoles sharing an edge with the seed point are the first-order neighborhood dipoles of the seed point, the dipoles sharing an edge with the first-order neighborhood dipoles and not sharing an edge with the seed point are the second-order neighborhood of the seed point, the dipoles sharing an edge with the m-order neighborhood dipoles and not sharing an edge with the m-1-order neighborhood dipoles are the m+1-order neighborhood dipoles of the seed point, 2≤m≤ N -1;

[0016] Steps S201 and S202 are repeated until the dipole regionalization is completed.

[0017] In some optional embodiments, the dipole intensity generalized total variation sparsity constraint model is wherein, is a penalty parameter, is a linear transformation matrix embodying the first-order variation characteristic of the source space dipole activation intensity, is a first-order variation value of the dipole intensity, , is a second-order variation value of the dipole intensity.

[0018] In some optional embodiments, the inverse problem solving model is wherein, is the dipole activation intensity to be solved, is an intermediate variable, which is an estimated first-order variation value of the dipole activation intensity, is a signal measured by a scalp electroencephalogram / magnetoencephalogram device, is the conduction matrix, is the number of electroencephalogram electrodes / magnetic sensors possessed by the scalp electroencephalogram / magnetoencephalogram device, , and is a penalty parameter, , , is a weight of the penalty term, , and is an infinitesimal quantity, , and are respectively , and are the mean values of

[0019] In some optional embodiments, the linear transformation matrix embodying the variation characteristic of the source space dipole activation intensity has the expression:

[0020]

[0021] wherein, , , , is the total number of edges of the dipoles in the source model.

[0022] In some optional embodiments, the determining the activation range and intensity of the dipoles in the source model based on the inverse problem solution model comprises:

[0023] deforming the inverse problem solution model into: , wherein, , and is a defined latent variable, is an auxiliary sparse matrix traversing all dipole groups in the source space, with a row number of and a column number of satisfying ;

[0024] introducing a scaled augmented Lagrange operator, and solving the deformed inverse problem solution model based on an alternating update manner to determine the activated dipoles and their intensity in the source model.

[0025] In some optional embodiments, the alternating update rule is:

[0026]

[0027] , are the serial numbers of the update times, is a Lagrange penalty operator, , and is the scaled augmented Lagrange operator, and is a soft threshold function of 1-norm, is a soft threshold function of 2-norm, and when is the independent variable, the above soft threshold function has the following mathematical expression:

[0028] , or ,

[0029] , .

[0030] In some optional embodiments, the end condition of the alternating update is: , is a first preset threshold, or, a first preset number threshold is reached.

[0031] In some optional embodiments, the solving the deformed inverse problem solving model based on the alternating updating comprises:

[0032] solving the deformed inverse problem solving model based on the updated , , , , , , , solving the deformed inverse problem solving model based on the updated

[0033] If the weight of the penalty term meets the preset condition or the current solving number reaches the second number threshold, the solving of the deformed inverse problem solving model is completed.

[0034] If the weight of the penalty term does not meet the preset condition and the current solving number has not reached the second number threshold, the , , , , , , , is updated again based on the updated , , , , , , , the deformed inverse problem solving model is solved again until the weight of the penalty term meets the preset condition or the solving number reaches the second number threshold based on the solving result.

[0035] In some optional embodiments, the preset condition required to be met by the weight of the penalty term is:

[0036] , , ;

[0037] wherein, THR1 is a second preset threshold, THR2 is a third preset threshold, THR3 is a fourth preset threshold, , , are respectively serial numbers of the updating numbers.

[0038] According to a second aspect, an embodiment of the present application provides a source imaging inverse problem analysis device, comprising:

[0039] a first model construction module, configured to construct a real brain model and a source model;

[0040] a calculation module, configured to solve a conduction matrix based on a structure of electrodes / sensors, the real brain model, the source model, and a positional relationship between the electrodes / sensors and the real brain model;

[0041] a second model construction module, configured to regionalize dipoles in the source model according to spatial distances, and construct an inter-group sparse constraint model;

[0042] a third model construction module, configured to construct a dipole intensity generalized total variation sparse constraint model;

[0043] a fourth model construction module, configured to fuse the inter-group sparse constraint model and the dipole intensity generalized total variation sparse constraint model, and construct an inverse problem solving model based on the conduction matrix;

[0044] a determination module, configured to determine a dipole activation range and intensity in the source model based on the inverse problem solving model.

[0045] According to a third aspect, an embodiment of the present application provides an electronic device, comprising:

[0046] a memory and a processor, which are in communication connection with each other, and the memory is configured to store a computer program, and the computer program is executed by the processor to implement any one of the source imaging inverse problem analysis methods according to the first aspect.

[0047] According to a fourth aspect, an embodiment of the present application provides a computer readable storage medium, configured to store a computer program, and the computer program is executed by a processor to implement any one of the source imaging inverse problem analysis methods according to the first aspect.

[0048] The embodiments of the present application provide a source imaging inverse problem analysis method and device based on group variation sparse constraint, and an electronic device. The embodiments of the present application first consider the cluster characteristics of neuron discharge, regionalize dipoles according to spatial distances of the dipoles, and consider inter-group sparsity to depict structural sparse prior of cortical electrical activity. In addition, the embodiments of the present application depict abruptness at an activation region boundary on the cortex and smoothness within the activation region through dipole intensity generalized total variation, to realize accurate estimation of dipole activation range and intensity. Finally, the embodiments of the present application realize accurate estimation of automatic activation range and intensity in a reweighted solving manner. BRIEF DESCRIPTION OF DRAWINGS

[0049] The features and advantages of the present application will be appreciated upon reference to the following detailed description and drawings, in which:

[0050] Figure 1 A flowchart of a source imaging inverse problem analysis method provided by an embodiment of the present application is shown in FIG. 2;

[0051] Figure 2 A process diagram of a source imaging inverse problem analysis method provided by an embodiment of the present application is shown in FIG. 3;

[0052] Figure 3 A comparison result diagram of a source imaging inverse problem analysis method provided by an embodiment of the present application and two existing methods is shown in FIG. 4;

[0053] Figure 4 A structure diagram of a source imaging inverse problem analysis device provided by an embodiment of the present application is shown in FIG. 5;

[0054] Figure 5 A structure diagram of an electronic device provided by an embodiment of the present application is shown in FIG. 6. DETAILED DESCRIPTION

[0055] In order to make the objectives, technical solutions and advantages of the embodiments of the present application clearer, the following will clearly and completely describe the technical solutions in the embodiments of the present application with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are some but not all of the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative efforts fall within the scope of the present application.

[0056] It should be noted that the terms “comprising”, “containing” or any other variant thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such a process, method, article or device. Without more limitations, the element defined by the phrase “comprising a” does not exclude the presence of additional identical elements in the process, method, article or device comprising the element. In addition, the terms “first”, “second” and the like are only for descriptive purposes and cannot be understood as indicating or implying relative importance or implicitly indicating the number of indicated technical features. In the description of the following embodiments, the meaning of “a plurality of” is two or more, unless otherwise specifically limited.

[0057] The goal of the inverse problem analysis method of source imaging is to accurately solve the activation range (position) and intensity of the cortical space (i.e. source space) dipole. The existing methods can be divided into methods based on the equivalent current dipole model and methods based on the distributed current dipole model. Among them, the method based on the equivalent current dipole model is to fix some parameters in the position, direction and intensity of the activated dipole under the condition of knowing the number of activated dipoles, adjust other unfixed parameters until the M / EEG generated by it matches the measured M / EEG. As can be seen, the equivalent current dipole model is a nonlinear model, and with the increase of the number of dipoles, its solution becomes complex. Generally, this method can only realize the position and intensity estimation of one or several dipoles. At the same time, it needs to take the number of dipoles as prior knowledge, which is difficult to accurately know in practice. In contrast, the distributed current dipole model is a linear model of cortical electrical activity and scalp EEG. It does not need to take the number of dipoles as prior knowledge, and is more in line with the actual situation.

[0058] The method based on the distributed current dipole model assumes that a large number of dipoles are uniformly distributed in the cortical space, and constructs the linear equation of the scalp electrode / sensor and the dipole according to the forward conduction model and solves the electrical activity in the cortical space. Since the number of scalp electrodes / sensors (≈10 2 ) is much smaller than the number of cortical space dipoles (≈10 4), source imaging is a highly underdetermined problem without a unique solution. Existing methods usually combine some constraints consistent with the characteristics of neuron firing to find an optimal solution. One norm and two norm are two important constraints in source imaging. Among them, two norm is a kind of constraint applied earlier, which takes into account the characteristics that neurons always work in an energy-saving way. This kind of method can obtain an analytical solution through the inverse operation of the matrix, and the running time is short. However, due to the constraint of the feasible region of two norm, for focal electrical activity, the two norm method usually obtains a very smooth solution, which makes the estimation of the activation range of the dipole inaccurate. Therefore, researchers introduce the one norm constraint which is more likely to obtain a sparse solution to source imaging, reducing the smoothness of the solution. A simplest method based on one norm constraint directly considers the sparsity of the number of activated dipoles, which can accurately estimate the center position of the activated dipole, but this constraint is too sparse to estimate the activation range of the dipole. Therefore, the edge sparsity is considered for solving, which uses the characteristics that the strength jump of adjacent dipoles only occurs on the boundary of the activated dipole and the non-activated dipole (first-order variation characteristics), to realize the estimation of the activation range of the activated dipole. However, this method only considers the sparsity of the edge, and does not consider the sparsity of the number of activated dipoles, when there are two clusters of activated dipoles with close distance, this method tends to solve a cluster of activated dipoles with a larger range. To solve the problem of the above method, the sparsity of the number of activated dipoles and the sparsity of the edge can be introduced at the same time, and a suitable penalty parameter is used to realize the estimation of the activation range of the dipole. However, the above edge sparsity constraint assumes that the dipole strength is constant in segments, which cannot accurately estimate the strength of the dipole in the activation range (i.e. the first-order variation characteristics have a step effect). In summary, an inverse problem analysis method of source imaging which can accurately solve the activation range and strength of the dipole in the cortex space is to be studied.

[0059] Please refer to Figure 1 and Figure 2 , the embodiment of the present application provides a kind of inverse problem analysis method of source imaging, comprising:

[0060] S101: construct real head model and source model;

[0061] Specifically, please refer to Figure 2A personalized real head model can be constructed based on fieldtrip according to the electrode / sensor structure of the scalp EEG device and the image obtained by magnetic resonance imaging (MRI); a source model can be constructed according to the MRI image by using software such as freesurfer; the source model can be a distributed source model, and the distributed source model means that the cortical discharge source (current dipole, referred to as 'dipole' for short) is uniformly distributed on the triangular meshed cortex; meanwhile, it is assumed that there is one dipole in a triangular mesh and the dipole is located at the center of gravity of the triangular mesh and the direction is the normal direction of the triangular mesh;

[0062] S102: a conductive matrix is solved based on the structure of the electrode / sensor, the real head model, the source model and the positional relationship between the electrode / sensor and the real head model;

[0063] S103: the dipoles in the source model are regionalized according to the spatial distance, and an inter-group sparse constraint model is constructed; here, based on the cluster discharge characteristics of neurons, it is considered that adjacent dipoles are more likely to be activated simultaneously, so the dipoles are regionalized according to the spatial distance.

[0064] S104: a dipole intensity generalized total variation sparse constraint model is constructed;

[0065] S105: the inter-group sparse constraint model and the dipole intensity generalized total variation sparse constraint model are fused, and an inverse problem solving model is constructed based on the conductive matrix;

[0066] S106: the activation range and intensity of the dipoles in the source model are determined based on the inverse problem solving model; specifically, the inverse problem solving model can be solved in a reweighted manner to automatically determine the activation range and intensity of the dipoles in the source model.

[0067] The embodiment of the application provides an inverse problem analysis method for source imaging based on group variation sparse constraint, which first considers the cluster characteristics of neuron discharge, regionalizes the dipoles according to the spatial distance of the dipoles, and considers the inter-group sparsity to depict the structural sparse prior of the cortical electrical activity. In addition, the dipole intensity generalized total variation is used to depict the abruptness at the boundary of the activated region on the cortex and the smoothness in the activated region, so as to realize accurate estimation of the activation range and intensity of the dipoles. Finally, the activation range and intensity are estimated automatically in a reweighted solving manner.

[0068] In some specific embodiments, the inter-group sparse constraint model is wherein, I is the total number of groups of the dipoles in the source model, is the number of dipoles in the i-th group, iA set of all dipole activation strengths in the group, a subset of a set of all elements in the source model, A set of all dipole activation strengths in the group, , The first, second, and so on, N s The activation strength of the first, second, and so on, N s The number of all dipoles in the source model, The index number of the group of dipoles in all dipoles, i The number of dipoles in the first The number of dipoles in the first i

[0069] In the embodiment of the application, when the dipoles in the source model are regionalized according to spatial distance, the dipoles in the source model are divided into I groups.

[0070] In some specific embodiments, the regionalization of the dipoles in the source model according to spatial distance comprises:

[0071] S201: Obtain the dipoles that have not been regionalized, and select one of the dipoles as a seed point;

[0072] S202: Select the dipoles in the n-order neighborhood of the seed point to form a group together with the seed point; wherein, , , The maximum order is The dipoles in the n-order neighborhood are all dipoles in the source model; the dipoles sharing an edge with the seed point are the first-order neighborhood dipoles of the seed point, the dipoles sharing an edge with the first-order neighborhood dipoles and not sharing an edge with the seed point are the second-order neighborhood dipoles of the seed point, the dipoles sharing an edge with the m-order neighborhood dipoles and not sharing an edge with the m-1-order neighborhood dipoles are the m+1-order neighborhood dipoles of the seed point, and 2≤m≤ N -1;

[0073] Repeat steps S201 and S202 until the regionalization of the dipoles is completed.

[0074] In the embodiment of the application, the spatial distance of a dipole is represented by the neighborhood of the dipole. When N is 10, the seed point and the first-order, second-order, third-order, fourth-order, fifth-order, sixth-order, seventh-order, eighth-order, ninth-order, and tenth-order neighborhood dipoles of the seed point form a dipole group.

[0075] ​In some embodiments, the dipole strength generalized total variation sparsity constraint model is wherein, is a penalty parameter, is a linear transformation matrix embodying the first order variation property of the source space dipole activation strength, is the first order variation value of the dipole strength, , is the second order variation value of the dipole strength.

[0076] Optionally, to ensure both the first order variation and the second order variation properties of the dipole strength, may be set to 0.5.

[0077] In some embodiments, the inverse problem solving model is wherein, is the dipole activation strength to be solved, is an intermediate variable, which is the estimated first order variation value of the dipole activation strength, is the signal measured by the scalp EEG / MEG device, is the conduction matrix, is the number of EEG electrodes / MEG sensors possessed by the scalp EEG / MEG device, , and is a penalty parameter, , , is the weight of the penalty term, , and is an infinitesimal quantity, , and are respectively , and are the mean values of ,

[0078] Optionally, , , , wherein, , , and are respectively the corresponding weights.

[0079] In some embodiments, the linear transformation matrix embodying the variation property of the source space dipole activation strength has the expression:

[0080]

[0081] wherein, , , , is the total number of edges of the dipoles in the source model.

[0082] In some specific embodiments, the determining the activation range and intensity of the dipoles in the source model based on the inverse problem solving model comprises:

[0083] deforming the inverse problem solving model into: , ; wherein, , and are defined latent variables, is an auxiliary sparse matrix traversing all dipole groups in the source space, the number of rows of the auxiliary sparse matrix is , the number of columns of the auxiliary sparse matrix is , and the auxiliary sparse matrix satisfies ;

[0084] introducing a scaled augmented Lagrangian operator, and solving the deformed inverse problem solving model based on an alternating update manner to determine the activated dipoles and the intensity thereof in the source model.

[0085] In the embodiments of the present application, the Alternating Direction Method of Multipliers (ADMM) method is used to solve the inverse problem solving model to obtain the activation of the dipoles.

[0086] In other alternative embodiments, a deep learning network can also be used to solve the above inverse problem solving model.

[0087] In some specific embodiments, the alternating update rule is:

[0088]

[0089] , are the serial numbers of the update times, is a Lagrangian penalty operator, , and are the scaled augmented Lagrangian operators; and are soft threshold functions of 1-norm, is a soft threshold function of 2-norm, and when is the independent variable, the above soft threshold function has the following mathematical expression:

[0090] , or ,

[0091] , .

[0092] In some alternative embodiments, the alternating update can be replaced by a deep neural network, such as a FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) network.

[0093] In some embodiments, the end condition of the alternating update is that: , is a first preset threshold, or a preset number threshold. The first preset threshold is a small number. That is, as long as one of the following conditions is met, the alternating update ends: is less than the first preset threshold, and the preset number threshold is reached.

[0094] For example, the end condition of the alternating update is that: or .

[0095] In some embodiments, the method of solving the deformed inverse problem solving model based on the alternating update comprises:

[0096] solving the deformed inverse problem solving model based on the updated , , , , , , , and determining whether the weight of the penalty term meets a preset condition and whether the current solving number reaches a second number threshold based on the solving result;

[0097] If the weight of the penalty term meets the preset condition or the current solving number has reached the second number threshold, the solving of the deformed inverse problem solving model is completed;

[0098] If the weight of the penalty term does not meet the preset condition and the current solving number has not reached the second number threshold, the , , , , , , , is updated again, and the deformed inverse problem solving model is solved based on the updated , , 、 、 、 、 、 Solving the deformed inverse problem solving model again until the weight of the penalty term is determined to meet the preset condition or the number of solving times reaches a second number threshold based on the solving result.

[0099] In some specific embodiments, the condition required to be met by the weight of the penalty term is:

[0100] 、 、 ;

[0101] Wherein, THR1 is a second preset threshold, THR2 is a third preset threshold, THR3 is a fourth preset threshold, , 、 are the serial numbers of the update times respectively.

[0102] In the embodiments of the present application, the weight of the penalty term is calculated based on the solving result of the inverse problem , and the solving of the inverse problem is repeated until the weight of the penalty term meets the preset condition or the number of solving times reaches a second number threshold.

[0103] For example, THR1, THR2 and THR3 can all be 10 -3 , and the second number threshold can be 10.

[0104] In this paper, a kind of classic inverse problem solving method related to the inverse problem solving method provided in the embodiments of the present application is selected, including: minimum one norm fusion variation method (TV- l 1, for example, source imaging based on structured sparsity method) and minimum norm difference fusion variation method (TV- l 1-2 , for example, sparsity and smoothness enhanced brain tomography method) as a comparative method. Figure 2 The inverse problem analysis method provided in the embodiments of the present application can also be called inverse problem solving method (Re-weighted Total Generalized Variation Fused Group Lasso, RTGVGS), and the positioning error diagram under the condition of determining dipole activation based on the same individualized head model and source model. Please refer to Figure 2 ​In this embodiment of the invention, a personalized mind model and a source model from a pre-established example of the fieldtrip public dataset are first obtained. Then, the dipoles in the source model are regionalized. Finally, an inverse problem-solving model with generalized total variational sparsity and group sparsity properties is constructed and solved.

[0105] To verify the accuracy of the proposed method, based on the obtained brain model source model, we constructed 50 simulated scalp EEGs with the same activation range under noise-free conditions and signal-to-noise ratios of 30dB, 20dB, 10dB, and 5dB, respectively. The simulated activation range was the 5th order neighborhood of the seed point dipole. The dipole activation intensity within the activation range was centered on the seed point dipole, with a mean of 0 and a standard deviation of [missing value]. ( The Gaussian distribution is the order of the dipole neighborhood of the seed point.

[0106] Dipole Localization Error (DLE) and Normalized Root Mean Square Error (NMSE) are used as performance evaluation metrics to assess the accuracy of the method in estimating the dipole activation range and intensity, respectively. DLE is calculated as follows:

[0107]

[0108] in, Represents DLE, , These represent the gold standard and the estimated number of activated dipoles, respectively. To activate the dipole as the gold standard, For the estimated activated dipole, For the first l The coordinates of the gold standard activated dipole For the first The estimated coordinates of the activated dipole.

[0109] The NMSE is calculated as follows:

[0110]

[0111] in, Represents NMSE, This represents the estimated dipole activation intensity. Based on... Figure 3It can be seen that, compared with the prior method, the source imaging inverse problem analysis method provided by the embodiment of the application can reduce the error of source space activated dipole positioning, more accurately estimate the activated strength of the dipole, and the result is less affected by noise. This shows that the generalized total variation sparsity and group sparsity prior in the model construction is more close to the characteristics of neuron discharge, and this method helps to estimate the activated state of the dipole, and also improves the robustness of the prior method.

[0112] Correspondingly, referring to Figure 4 The embodiment of the application provides a source imaging inverse problem analysis device, which comprises:

[0113] A first model construction module 301 is configured to construct a real brain model and a source model;

[0114] A calculation module 302 is configured to solve a conduction matrix based on the structure of electrodes / sensors, the real brain model, the source model, and the positional relationship between the electrodes / sensors and the real brain model.

[0115] A second model construction module 303 is configured to regionalize the dipoles in the source model according to spatial distance, and construct a group sparsity constraint model;

[0116] A third model construction module 304 is configured to construct a dipole intensity generalized total variation sparsity constraint model;

[0117] A fourth model construction module 305 is configured to fuse the group sparsity constraint model and the dipole intensity generalized total variation sparsity constraint model, and construct an inverse problem solving model based on the conduction matrix;

[0118] A determination module 306 is configured to determine the activated range and intensity of the dipoles in the source model based on the inverse problem solving model.

[0119] The embodiment of the application provides a source imaging inverse problem analysis device based on group variation sparsity constraint. The device first considers the cluster characteristics of neuron discharge, regionalizes the dipoles according to the spatial distance of the dipoles, and considers the group sparsity to depict the structural sparsity prior of the cortical electrical activity. In addition, the generalized total variation of the dipole intensity is used to depict the abruptness at the boundary of the activated region on the cortex and the smoothness in the activated region, so as to accurately estimate the activated range and intensity of the dipole. Finally, the automatic activated range and intensity estimation is realized by using the reweighted solving method.

[0120] In some specific embodiments, the group sparsity constraint model is wherein, I is the total number of groups of the dipoles in the source model, is the first group of dipoles in the source model, iThe set of all dipole activation intensities in the group is the vector of all dipole intensities in the source model. A subset of the set consisting of all elements in. , They are respectively the 1st, 2nd, ... N s The activation intensity of each dipole N s The number of all dipoles in the source model. For the first i A vector formed by the index numbers of the group of dipoles among all dipoles. For the first i The number of dipoles in the group.

[0121] In some specific implementations, the second model construction module 303 includes:

[0122] The seed point selection unit is used to obtain the dipoles that are currently not regionalized, and select one of the dipoles as a seed point.

[0123] A grouping unit is used to select dipoles in the nth-order neighborhood of the seed point and form a group with the seed point; wherein... , , It is the largest order, that is The dipoles within the first-order neighborhood are all the dipoles in the source model; the dipoles sharing an edge with the seed point are the first-order neighborhood dipoles of the seed point; the dipoles sharing an edge with the first-order neighborhood dipoles but not with the seed point are the second-order neighborhoods of the seed point; the dipoles sharing an edge with the m-order neighborhood dipoles but not with the m-1-order neighborhood dipoles are the m+1-order neighborhood dipoles of the seed point, where 2 ≤ m ≤ N -1;

[0124] The control unit is used to control the seed point selection unit and the grouping unit to execute repeatedly until the dipole regionalization is completed.

[0125] In some specific implementations, the generalized total variational sparse constraint model for dipole strength is as follows: ,in, For penalty parameters, To represent the linear transformation matrix with first-order variational properties of the activation intensity of the source space dipole, The first variational value of the dipole intensity. , This represents the second variational value of the dipole intensity.

[0126] In some specific implementations, the inverse problem solution model is as follows: ,in, is the dipole activation strength to be solved, is an intermediate variable, which is the first order variational value of the estimated dipole activation strength, is the signal measured by the scalp EEG / MEG device, is the conduction matrix, is the number of EEG electrodes / MEG sensors possessed by the scalp EEG / MEG device, , and is a penalty parameter, , , is the weight of the penalty term, , and is an infinitesimal quantity, , and are respectively , and the mean of ,

[0127] In some specific embodiments, the linear transformation matrix which embodies the variational characteristics of the source space dipole activation strength has the expression:

[0128]

[0129] wherein, , , , is the total edge number of the dipoles in the source model.

[0130] In some specific embodiments, the determining module 306 comprises:

[0131] a deformation unit, configured to deform the inverse problem solving model into:

[0132] , ; wherein, , and are defined latent variables, is an auxiliary sparse matrix traversing all dipole groups in the source space, with the number of rows being and the number of columns being , satisfying ;

[0133] a solving unit, configured to introduce the scaled augmented Lagrange operator, and solve the deformed inverse problem solving model based on an alternating updating manner to determine the activated dipoles and their strengths in the source model.

[0134] In some specific embodiments, the alternately updated rule is:

[0135]

[0136] , is the serial number of the update times, is a Lagrange punishment operator, , and is the scaled augmented Lagrange operator, and is a soft threshold function of 1-norm, is a soft threshold function of 2-norm, and is the independent variable, the above soft threshold function has the following mathematical expression:

[0137] , or ,

[0138] , .

[0139] In some optional specific embodiments, the end condition of the alternately updating is: , is a first preset threshold, or a first preset number threshold is reached.

[0140] In some specific embodiments, the solving unit is specifically configured to:

[0141] based on the updated , , , , , , , solve the deformed inverse problem solving model, and determine whether the weight of the penalty term meets a preset condition and whether the current solving number reaches a second number threshold based on a solving result;

[0142] if the weight of the penalty term meets the preset condition or the current solving number has reached the second number threshold, the solving of the deformed inverse problem solving model is completed;

[0143] if the weight of the penalty term does not meet the preset condition and the current solving number has not reached the second number threshold, the , is updated again. 、 、 、 、 、 , and based on the updated 、 、 、 、 、 、 、 solving the deformed inverse problem solving model again until it is judged based on the solving result that the weight of the penalty term meets a preset condition or the number of solving reaches a second number threshold.

[0144] In some specific embodiments, the condition required to be met by the weight of the penalty term is:

[0145] 、 、 ;

[0146] wherein THR1 is a second preset threshold, THR2 is a third preset threshold, THR3 is a fourth preset threshold, , 、 are respectively serial numbers of the number of updating.

[0147] The embodiment of the application is a device embodiment based on the same inventive concept as the above-mentioned method embodiment, and therefore specific technical details and corresponding technical effects are referred to the above-mentioned method embodiment, which will not be described here again.

[0148] The embodiment of the application further provides an electronic device, as shown in Figure 5 may include a processor 41 and a memory 42, wherein the processor 41 and the memory 42 can be connected to each other in communication through a bus or other means, Figure 5 for example, through a bus connection.

[0149] The processor 41 can be a central processing unit (CPU). The processor 41 can also be other general-purpose processors, digital signal processors (DSP), application specific integrated circuits (ASIC), field programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc. chips, or combinations of the above various types of chips.

[0150] Memory 42, as a non-transitory computer-readable storage medium, can be used to store non-transitory software programs, non-transitory computer-executable programs, and modules, such as the program instructions / modules corresponding to the source imaging inverse problem analysis method in this embodiment of the invention (e.g., Figure 4 The first model building module 301, the calculation module 302, the second model building module 303, the third model building module 304, the fourth model building module 305, and the determination module 306 are shown. The processor 41 executes various functional applications and data processing by running non-transitory software programs, instructions, and modules stored in the memory 42, thereby realizing the inverse problem analysis method of source imaging in the above method embodiment.

[0151] The memory 42 may include a program storage area and a data storage area. The program storage area may store the operating system and applications required for at least one function; the data storage area may store data created by the processor 41, etc. Furthermore, the memory 42 may include high-speed random access memory and may also include non-transitory memory, such as at least one disk storage device, flash memory device, or other non-transitory solid-state storage device. In some embodiments, the memory 42 may optionally include memory remotely located relative to the processor 41, and these remote memories may be connected to the processor 41 via a network. Examples of such networks include, but are not limited to, the Internet, corporate intranets, local area networks, mobile communication networks, and combinations thereof.

[0152] The one or more modules are stored in the memory 42, and when executed by the processor 41, they perform actions such as... Figures 1-2 The inverse problem analysis method for source imaging in the illustrated embodiment.

[0153] For specific details regarding the aforementioned electronic devices, please refer to the relevant documentation. Figures 1 to 2 The relevant descriptions and effects in the illustrated embodiments are for understanding purposes only and will not be repeated here.

[0154] Accordingly, this embodiment of the invention also provides a computer-readable storage medium for storing a computer program. When the computer program is executed by a processor, it implements the various processes of the above-described source imaging inverse problem analysis method embodiment and achieves the same technical effect. To avoid repetition, it will not be described again here.

[0155] Computer-readable media includes permanent and non-permanent, movable and non-movable media that can be implemented by any method or technology to store information. The information can be computer-readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, compact disc read-only memory (CD-ROM), digital versatile disc (DVD) or other optical storage, magnetic cassette, magnetic tape disk storage or other magnetic storage devices, or any other non-transmission medium that can be used to store information accessible by a computing device. According to the definition herein, computer-readable media does not include transitory media such as modulated data signals and carriers.

[0156] Each of the embodiments in the specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the difference from other embodiments. In particular, for the device embodiments, since they are basically similar to the method embodiments, the description is relatively simple, and the relevant parts can be referred to the part of the method embodiments.

[0157] The above only describes the embodiments of the present application and does not limit the present application. Those skilled in the art can make various changes and modifications to the present application. Any modification, equivalent replacement, improvement, etc. within the spirit and principle of the present application shall be included in the scope of the claims of the present application.

Claims

1. A method of source imaging inverse problem analysis, characterized by, The method comprises the following steps: constructing a real brain model and a source model; solving a conduction matrix based on the structure of electrodes / sensors, the real brain model, the source model, and the position relationship between the electrodes / sensors and the real brain model; regionalizing dipoles in the source model according to spatial distances, and constructing an inter-group sparse constraint model; constructing a dipole intensity generalized total variation sparse constraint model; fusing the inter-group sparse constraint model and the dipole intensity generalized total variation sparse constraint model, and constructing an inverse problem solving model based on the conduction matrix; determining the activation range and intensity of dipoles in the source model based on the inverse problem solving model; Wherein, the dipole intensity generalized total variation sparse constraint model is Wherein, is a penalty parameter, is a linear transformation matrix embodying that the source space dipole activation intensity has a first-order variation characteristic, is a first-order variation value of the dipole intensity, , is a second-order variation value of the dipole intensity.

2. The method of claim 1, wherein, The inter-group sparse constraint model is wherein, I is the total number of groups of dipoles in the source model, is the set of activation intensities of all dipoles in the i th group, is the subset of the set of all elements in the vector of activation intensities of all dipoles in the source model, , is the activation intensity of the 1st, 2nd, …, th dipole, respectively, N s N s is the total number of dipoles in the source model, is the vector of index numbers of the dipoles in the i th group among all dipoles, is the number of dipoles in the i th group.​ 3. The method of claim 1, wherein, the regionalizing dipoles in the source model according to spatial distances comprises: S201: obtaining the dipoles that have not been regionalized, and selecting one of the dipoles as a seed point; S202: selecting dipoles in the n-order neighborhood of the seed point to form a group together with the seed point; wherein, , , is the maximum order, that is, dipoles in the n-order neighborhood are all dipoles in the source model; dipoles sharing an edge with the seed point are first-order neighborhood dipoles of the seed point, dipoles sharing an edge with the first-order neighborhood dipoles and not sharing an edge with the seed point are second-order neighborhood dipoles of the seed point, dipoles sharing an edge with the m-order neighborhood dipoles and not sharing an edge with the m-1-order neighborhood dipoles are (m+1)-order neighborhood dipoles of the seed point, 2≤m≤ N -1; repeating steps S201 and S202 until the regionalization of the dipoles is completed.

4. The method of claim 2, wherein, The inverse problem solving model is wherein, is the dipole activation strength to be solved, is an intermediate variable, which is the first order variation value of the estimated dipole activation strength, is the signal measured by the scalp electroencephalogram / magnetoencephalogram device, is the conduction matrix, is the number of electroencephalogram electrodes / magnetic sensors possessed by the scalp electroencephalogram / magnetoencephalogram device, , and is a penalty parameter, , , is the weight of the penalty term, , and is an infinitesimal quantity, , and are the mean values of , and respectively.

5. The method of claim 2, wherein, The linear transformation matrix embodying the source spatial dipole activation strength has a variational characteristic The expression is: wherein , , , is the total number of edges of the dipoles in the source model.

6. The method of claim 4, wherein, the determining the activation range and intensity of dipoles in the source model based on the inverse problem solving model comprises: The inverse problem solving model is deformed as: , ; wherein, , and are defined latent variables, is an auxiliary sparse matrix traversing all dipole groups in the source space, with rows and columns, satisfying ; introducing a scaled augmented Lagrange operator, and solving the deformed inverse problem solving model based on an alternating updating manner to determine the activated dipoles and their intensities in the source model.

7. The method of claim 6, wherein, the alternating updating rule is: , These are the sequence numbers of the update counts. To punish operators for Lagrange, , and For the scaled augmented Lagrange operator, and A soft thresholding function with a 1-norm. The soft thresholding function is a 2-norm function. When the variable is 'soft threshold', the mathematical expression for the soft threshold function is as follows: , or , ， 。 8. The method according to claim 6 or 7, characterized in that, The ending condition of the alternately updating is: , is the first preset threshold, or, a first preset number threshold is reached.

9. The method of claim 7, wherein, the solving the deformed inverse problem solving model based on the alternating updating manner comprises: based on the updated , , , , , , , determines whether the weight of the penalty term satisfies a preset condition and whether the current number of solutions reaches a second number threshold. if the weight of the penalty term meets the preset condition or the current solving number has reached the second number threshold, the solving of the deformed inverse problem solving model is completed; If the weight of the penalty term does not satisfy the preset condition and the current number of solving has not reached the second number threshold, the weight of the penalty term is updated again 、 、 、 、 、 、 、 , and based on the updated 、 、 、 、 、 、 、 The weight of the penalty term is calculated again until the preset condition is satisfied or the number of solving reaches the second number threshold.

10. The method of claim 9, wherein, the preset condition required to be met by the weight of the penalty term is: 、 、 ； Wherein, THR1 is the second preset threshold, THR2 is the third preset threshold, THR3 is the fourth preset threshold, are respectively serial numbers of the update times.​​ 11. A source imaging inverse problem analysis apparatus characterized by comprising: The method comprises the following steps: a first model construction module for constructing a real brain model and a source model; a calculation module for solving a conduction matrix based on the structure of electrodes / sensors, the real brain model, the source model, and the position relationship between the electrodes / sensors and the real brain model; a second model construction module for regionalizing dipoles in the source model according to spatial distances, and constructing an inter-group sparse constraint model; a third model construction module for constructing a dipole intensity generalized total variation sparse constraint model; a fourth model construction module for fusing the inter-group sparse constraint model and the dipole intensity generalized total variation sparse constraint model, and constructing an inverse problem solving model based on the conduction matrix; a determination module for determining the activation range and intensity of dipoles in the source model based on the inverse problem solving model; Wherein, the dipole intensity generalized total variation sparse constraint model is Wherein, is a penalty parameter, is a linear transformation matrix embodying that the source space dipole activation intensity has a first-order variation characteristic, is a first-order variation value of the dipole intensity, , is a second-order variation value of the dipole intensity.

12. An electronic device, comprising: The method comprises the following steps: a memory and a processor, which are communicatively connected with each other, the memory is used for storing a computer program, and the computer program is executed by the processor to implement the inverse problem analysis method for source imaging in any one of claims 1 to 10.

13. A computer-readable storage medium, characterized in that, The computer readable storage medium is used for storing a computer program, and the computer program is executed by the processor to implement the inverse problem analysis method for source imaging in any one of claims 1 to 10.

Images