CN1903119A - Electrocerebral source positioning method based on discrete restraint - Google Patents

Abstract

An electrocerebral source locating method based on sparse constraint includes such steps as determining transfer array A, obtaining actually recorded electrocerebral signal Y via multi-channel electrocerebral signal recording system, pre-processing to determine the time when the source is analyzed, initializing source vector, iteration stop error and maximal iteration times, updating the diagonal array of weight, applying sparse constraint, updating source information vector, judging the condition for stopping iteration, and stopping or continuous iteration.

Description

Brain power supply localization method based on the constraint of sparse property

Technical field:

The invention belongs to the biology information technology field, relate to a kind of method that the brain power supply is positioned, be mainly used in human brain function and with the research and the diagnosis of human brain relevant disease.

Background technology:

According to the space operation position of apparent survey current potential inverting location brain power supply is a major issue in the EEG research, and it is a nonlinear optimization inverse problem in essence.In order to simplify complexity of the calculation, in the inverting location of brain power supply, a linear method commonly used goes to approach nonlinear problem, thereby can be summed up as following linear system (C.M.Michel, a M.Murray to electroencephalography (eeg) inverse problem, EEG source imaging, Clinical neurophysiology, 115,2195-2222,1997)

b＝Ax (1)

Wherein b is the vector of a M * 1 for the current potential that head table electrode records, and M is a table record electrode number; A is transfer matrix (Lead Field), is the matrix of one M * N, and N is the dimension of source activity solution space; X is the vector of a N * 1, is pending sterically defined source information vector.In electroencephalography (eeg) inverse problem, common M＜＜N, therefore, the linear system shown in the formula (1) is one and seriously owes fixed system.To top equation, there are infinite many groups to separate x and satisfy the voltage distribution b that observes, so in order to obtain separating of realistic physiological property, people often according to different needs and research purpose to top inverse problem apply suitable constraint (Yao is German-Chinese. electricity theory and method that brain function is surveyed. Beijing: Science Press, 2003,195-243).In the research in early days, mainly adopted minimum modulus to separate (minimum norm solution, MNS), develop into current weighting mould commonly used on this basis and separated (weightedminimum-norm solution, WMN), wherein the most representative is LORETA (Lowresolution electromanetic tomography) method, and it is the method that a kind of Laplacian of employing differential operator is weighted solution space.Adopt the WMN method, usually can only obtain the active block fuzzy region in source, do not reach the meticulousr source location (R.D.Pascual-Marqui that Neuroscience Research needs, Review of methods for solving the EEG inverseproblem, Int.J Bioelectromagnetism, 1,75-86,1999; C.M.Michel, M.Murray, EEG source imaging, Clinical neurophysiology, 115,2195-2222,1997)

Because human brain is when carrying out certain task or having certain disease, the main neural activity zone of its correspondence is locality and sparse property, therefore, if the constraint when these two characteristics are found the solution as under determined system can obtain to meet more (the C.M.Michel that separates of brain power supply ambulatory physiological characteristic in theory, M.Murray, EEG source imaging, Clinicalneurophysiology, 115,2195-2222,1997; Yao is German-Chinese. electricity theory and method that brain function is surveyed. and Beijing: Science Press, 2003,195-243).In current electroencephalography (eeg) inverse problem technology, set up the sparse solution that dual mode obtains the source: the one, be called as a kind of iteration WMN method of FOCUSS (Focal underdetermined system solver) algorithm; The 2nd, directly find the solution l _pSeparating of (p≤1) modular constraint.FOCUSS passes through progressively iteration, make most of element trends zero of solution space, thereby make the energy localization (I.F.Gorodnitsky that separates, B.D.Rao, Sparse signal reconstruction from limited datausing FOCUSS:A re-weighted minimum norm algorithm, IEEE Trans.S.P., 45,600-616,1997).The further iteration of fuzzy solution that it normally obtains methods such as LORETA is obtained the sparse solution that energy localizes.The primary iteration process of FOCUSS can be stated as: adopt a linear transformations x=Wq, the weighting matrix W of current iteration _kBe taken as by a preceding iteration result element x _K-1The diagonal matrix that constitutes is designated as W _k=(diag (x _K-1)), q is an intermediate variable.FOCUSS can be finished by following 3 step iteration:

1 ⁰ W _k＝(diag(x _k-1))

2 ⁰ q _k＝(AW _k) ⁺b

3 ⁰ x _k＝W _kq _k

(AW wherein _k) ⁺Expression is to matrix A W _kCarry out generalized inverse finding the solution.

From the iterative process here as can be seen, a process of finding the solution matrix inversion is arranged in the iterative process of FOCUSS,, determined the stability of FOCUSS algorithm to a great extent and the immunocompetence of noise to the estimation quality of matrix inversion.Current improvement to FOCUSS mainly is with various regularization technology, comprises that singular value blocks the characteristic of improving matrix inversion with various forms of Tikhonov regularization technology.The 2nd kind is adopted l _pThe method of (p≤1) modular constraint, the inverse problem that is exactly wushu (1) expression is converted into following problem (C.Silva, J.C.Maltez, E.Trindade, Evaluation of L1 and L2 minimum norm performances on EEGlocalizations, Clinical neurophysiology, 115,1657-1668,2004)

arg min‖b-Ax‖ ₂+λ‖x‖ _p (2)

Directly adopt optimization method to find the solution then, to obtain a sparse source distribution x to (2) formula.In electroencephalography (eeg) inverse problem, common p=1.Here, the constraint of sparse property is to introduce by an item that is similar to the Tikhonov regularization, and the effect of its existing sparse property constraint also has the effect of regularization, thereby makes that whole computational process is sane.But this technology does not have iterative process, and the sparse property of separating that it is obtained is not so good as the FOCUSS method.

Above-mentioned FOCUSS and l _p(p≤1) modulo n arithmetic all is the sparse solution in the source of asking for, but two kinds of technology respectively have characteristics.Present technique is based oneself upon the active sparse property physiological property of brain power supply, with l _p(p≤1) modular constraint is fused among the iterative process of FOCUSS, has obtained more sane sparse source location result.

Summary of the invention

Technical problem to be solved by this invention is to provide a kind of based on l _pThe iteration WMN brain power supply localization method of modular constraint (being called for short LIPSS).This method is with l _pThe constraint of the sparse property of mould is fused among the iterative process of FOCUSS, to obtain more stable and sparse source location result.

The technical scheme that the present invention solve the technical problem employing is that the brain power supply localization method based on the constraint of sparse property may further comprise the steps:

1) determines transfer matrix A;

2) obtain the EEG signals Y of physical record by multiple tracks EEG signals recording system, carry out pretreatment, determine to carry out the moment that the source is analyzed;

3) initialization source vector x _K-1, k=1, iteration stops the iterations of an error ε and a maximum;

4) upgrade weight diagonal matrix: W _k=(diag (x _K-1));

5) add the constraint of sparse property, utilize optimization method direct estimation q _k: arg min ‖ b-AW _kq _k‖ ₂+ λ ‖ q _k‖ _p

6) upgrade source information vector: x _k=W _kq _k

7) stopping criterion for iteration is judged: the source distribution before and after relatively upgrading changes, as ‖ x _k-x _K-1When ‖≤ε or iterations exceeded the maximum iteration time restriction of setting, iteration stopped, x _kBe the final location estimation result in source; Otherwise k=k+1 changes step 2), continue iteration;

Wherein, λ ≠ 0, W _k≠ I (I is a unit matrix).

In the described step 1, comprise step by step following:

A1), the head of object to be measured is carried out MRI or CT scan, obtain the image information of cranial anatomy structure;

A2), extraction step a1) brain part in the image information of gained, then brain is cut apart, extract the functional areas, source (mainly comprising positions such as grey matter, Hippocampus, cerebellum) of brain part again;

A3), with the grid of certain precision with step a2) functional areas, brain source of gained carry out subdivision, determine solution space grid (the locus sequence number that comprises dimension He each grid of solution space);

A4), determine the spatial positional information of each electrode of multiple tracks EEG signals recording system;

A5), determine the model of brain power supply;

A6), utilize step a3) to step a5) in solution space grid, electrode position information and the brain power source model determined, utilize forward modeling method to calculate transfer matrix A, concrete grammar is as follows: the source of placing unit on each solution space position, utilize numerical computation method to calculate the Potential distribution that this unit source produces at the electrode position place, this Potential distribution constitutes the string in the transfer matrix, by that analogy, after all solution space traversals are placed the unit source, just can obtain transfer matrix A.

Especially, the p value is 1.

Described step 2) pretreatment in comprises: filtering, baseline correction, eye electricity are removed.

The invention has the beneficial effects as follows, the former method of comparing, the present invention adopts the LPISS technology, from the active sparse property physiology fact of brain power supply, with l _p(p≤1) modular constraint is fused in the iterative process of FOCUSS, can obtain more stable sparse source location result.

The present invention is further illustrated below in conjunction with the drawings and specific embodiments.

Description of drawings

Fig. 1: under NSR=15% to LORETA, FOCUSS and the LPISS positioning result sketch map in two sources.Chromatic square region is being represented the source position of estimating, different colors is being represented different energy; Sign has the colored square region of blue spider representing to estimate lap between source position and dummy source position; Blue spider in green circle is represented the not lap between dummy source and estimation source.

Wherein a is the LORETA positioning result, a left side (23.00 ,-51.00,51.67) (mm), right (43.00 ,-41.00,61.67) (mm).

B is the FOCUSS positioning result, a left side (23.00 ,-51.00,51.67) (mm), right (43.00 ,-41.00,61.67) (mm).

C is the LPISS positioning result, a left side (23.00 ,-51.00,51.67) (mm), right (43.00 ,-41.00,61.67) (mm).

Fig. 2: under the different noise conditions to the positioning performance indicatrix in two sources.Wherein a is source 1, and b is source 2.

Fig. 3: the depth is mixed the positioning result sketch map in source.Chromatic square region is being represented the source position of estimating, different colors is being represented different energy; Sign has the colored square region of blue spider representing to estimate lap between source position and dummy source position; Blue spider in green circle is being represented the not lap between dummy source and estimation source.Wherein a is the LORETA positioning result, a1:(-83.00,19.00,31.67) (mm); A2:(-23.00,29.00 ,-8.33) (mm); A3:(-63.00,39.00 ,-48.33) (mm).

B is the FOCUSS positioning result, b1:(-83.00,19.00,31.67) (mm); B2:(-23.00,29.00 ,-8.33) (mm); B3:(-63.00,39.00 ,-48.33) (mm).

C is the LPISS positioning result.c1：(-83.00，19.00，31.67)(mm)；c2：(-23.00，29.00，-8.33)(mm)；c3：(-63.00，39.00，-48.33)(mm)。

Fig. 4: LPISS is to the source location result schematic diagram of IOR data.Wherein a is right frontal lobe eye district, and b is a thalamus, and c is right cerebellum, and d is a left occipital lobe.

The specific embodiment

The present invention adopts the LPISS technology, specifically, may further comprise the steps:

Step 1. is determined transfer matrix A, comprises step by step following:

A1) head to object to be measured carries out MRI or CT scan, obtains the image information of cranial anatomy structure;

A2) extraction step a1) the brain part in the image information of gained is cut apart brain then, extracts the functional areas, source (mainly comprising positions such as grey matter, Hippocampus, cerebellum) of brain part again;

A3) with the grid of certain precision with step a2) functional areas, brain source of gained carry out subdivision, determine solution space grid (the locus sequence number that comprises dimension He each grid of solution space);

A4) determine the spatial positional information of each electrode of multiple tracks EEG signals recording system;

A5) determine the model of brain power supply;

A6) utilize step a3) to step a5) in solution space grid, electrode position information and the brain power source model determined, utilize forward modeling method to calculate transfer matrix A, concrete grammar is as follows: the source of placing unit on each solution space position, utilize numerical computation method to calculate the Potential distribution that this unit source produces at the electrode position place, this Potential distribution constitutes the string in the transfer matrix, by that analogy, after all solution space traversals are placed the unit source, just can obtain transfer matrix A;

Step 2. is obtained the EEG signals Y of physical record by multiple tracks EEG signals recording system, and carries out necessary as basic pretreatment such as filtering, baseline correction, removals of eye electricity, and the moment of determining to carry out the source analysis;

Step 3. initialization source vector x _K-1, k=1, iteration stops the iterations of an error ε and a maximum;

Step 4. is upgraded weight diagonal matrix: W _k=(diag (x _K-1));

Step 5. adds the constraint of sparse property, utilizes optimization method direct estimation q _k:

arg min‖b-AW _kq _k‖ ₂+λ‖q _k‖ _p；

Step 6. is upgraded source information vector: x _k=W _kq _k

Step 7. stopping criterion for iteration is judged: the source distribution before and after relatively upgrading changes, as ‖ x _k-x _K-1When ‖≤ε or iterations exceeded the maximum iteration time restriction of setting, iteration stopped, x _kBe the final location estimation result in source; Otherwise k=k+1 changes step 2, continues iteration.

Preferred value is p=1; In the 3rd step, adopt and become the optimization method of yardstick l based on BFGS _p(p≤1) modulus problem is found the solution.

The grid of the certain precision in the such scheme, the step a3 of step 1.) is taken all factors into consideration computational accuracy and efficient, generally gets the 10mm/ lattice; The step a4 of step 1.) the multiple tracks EEG signals recording system described in can be the EEG signals recording system of 32 roads, 64 roads, 128 roads and the 256 road electrodes of standard; The step a5 of step 1.) the brain power source model described in is generally point charge model or dipole model; The step a6 of step 1.) numerical computation method described in can be boundary element algorithm or finite element algorithm; The initial solution in the source in the step 3 is made as separating of LORETA; P=1 in the step 5 adopts based on BFGS to become the optimization method of yardstick to l _p(p≤1) modulus problem is found the solution.Step 1 in fact can be according to existing techniques in realizing.

As effect, based on a real standard header model, utilize dummy source that algorithm has been carried out the positioning performance analysis relatively, and be applied to returning inhibition (Inhibition of return, IOR) the source location analysis of brain electricity experimental data.Transfer matrix A specifically calculates in such a way: by the MRI head model of scanning acquisition, the dipole source moving position is limited to grey matter, Hippocampus and other possible zones of action, source of brain, be separated into 910 positions by the 10mm mesh generation, on each grid position, place respectively along X, Y, the dipole of Z direction, employing standard 128 road electrode systems, (Boundary Elements Method, BEM) method is calculated and is obtained transfer matrix A to utilize boundary element.When implementing, the maximum iteration time of LPISS is 100, and it is 1E-5 that algorithm iteration stops error ε.

When the assessment algorithm performance, three following performance indications have mainly been adopted: 1) space orientation error, E _Localization, be defined as dummy source and in surround a bulbous region of dummy source, have the distance between the position in estimation source of ceiling capacity; 2) source energy error, E _Energy, can be used for assessment algorithm to the recovery capability of source energy, its calculating is defined as: $E_{\mathrm{energy}}=|1-\frac{\left|\right|J_{\max }\left|\right|}{\left|\right|J_{\mathrm{simu}}\left|\right|}|,$ J wherein _SimuBe the dipole moment of dummy source, J _MaxIt is the dipole moment (C.M.Michel, M.Murray, EEG source imaging, Clinical neurophysiology, 115,2195-2222,1997) of the dipole in the spheric neighbo(u)rhood of corresponding dummy source certain radius on every side with ceiling capacity; 3) at area-of-interest (region of interest, ROI) fuzzy index (the normalized burring index of Nei normalization, NBI), NBI can be used for the spatial resolving power of metric algorithm, NBI is defined as (D Yao, B.He, A Self-CoherenceEnhancement Algorithm and its Application to Enhancing 3D SourceEstimation from EEGs, Annals of Biomedical Engineering, 29 (11), 1019-1027,2001)

${\mathrm{NBI}}_i=\sqrt{\underset{k}{\mathrm{\Σ}}{\left|\right|r_k-r_i\left|\right|}^2{J}^2\left(k\right)}/\sqrt{\underset{k}{\mathrm{\Σ}}{\left|\right|r_k-r_i\left|\right|}^2/\underset{k}{\mathrm{\Σ}}1}$

In this work, when calculating NBI, employing all be spherical ROI.Each symbolic significance in NBI definition is as follows: subscript i is representing the grid sequence number of central point of the spherical ROI of NBI to be calculated, and when the NBI of the source distribution of simulation hypothesis was calculated, i was exactly the pairing grid sequence number of true spatial location of dummy source; When the NBI that estimates source distribution was calculated, i elected the pairing grid sequence number of the estimation dipole position that has ceiling capacity in the certain spherical ROI of dummy source as.Subscript k is meant all the effective grid node ID in ROI.r _k, r _iBe respectively that sequence number is k, the pairing locus of the mesh point of i coordinate.NBI can be characterized in the distribution situation in the source in the ROI zone.Usually, the source distributes mild more in ROI, and NBI approaches 1 more; If present sharp-pointed distribution in that ROI is endogenous, then NBI approaches 0.Obviously, the source of Discrete Distribution has less NBI.An outstanding source location algorithm should have less E _LocalizationAnd E _Energy, the NBI of the source distribution of Gu Jiing should be more similar with the real source NBI that distributes simultaneously.

1. the simulation experiment of positioning performance is analyzed relatively:

Method: at first tested the source location ability of algorithm under different noise conditions.Two dipole moments are respectively the dipole of (1.90,0.60,0.40) and (0.40,2.30,0.00), be placed on two insular position (73.00 ,-51.00,31.67) (mm) and (43.00 ,-41.00,61.67) (mm).Apparent survey current potential is just to drill calculating by BEM, and obtains after applying the Gaussian noise of different N SR.To different NSR, adopt LORETA respectively, FOCUSS and LPISS position the source.In this work, NSR is defined as the ratio of noise criteria difference and signal standards difference.Calculate E _Energy, E _LocalizationThe radius of the ROI that adopts during with NBI all is 25mm, and different is to calculate E _Energy, E _LocalizationThe time the center of spherical ROI in the dummy source position, the center of the spherical ROI when calculating NBI is (dipole of ceiling capacity) position in the estimation source.Positioning result when NSR=15% is presented in the accompanying drawing 1, and accompanying drawing 2 has shown the positioning performance index under all analogue noise situations.

Next has been tested algorithm and the depth has been mixed the positioning result in source.In this simulation experiment, the dipole moment intensity of three dipoles is respectively (6.00,2.70,1.90), (8.00,3.00,1.00) and (5.80,1.00,2.00), they are placed on three discrete positions (83.00 respectively, 19.00,31.67) and (mm), (23.00,29.00 ,-8.33) and (mm) with (63.00,39.00,-48.33) (mm), wherein first is two shallow sources with the 3rd source, and second is deep focus.Equally, also be to obtain head table Potential distribution, and apply the Gaussian noise of NSR=15% to it with the BEM forward modeling method.Also be to utilize LORETA, FOCUSS and three kinds of methods of LPISS that it is positioned, calculate corresponding separately E by top same method _Localization, E _EnergyWith the NBI index.The result is presented in accompanying drawing 3 and the table 1.

Table 1. three source location results' under NSR=15% performance indications

	Source 1	Source 2	Source 3
				LORETA FOCUSS LPISS	LORETA FOCUSS LPISS	LORETA FOCUSS LPISS
	E energy	0.8980 0.4419 0.2502	0.9749 0.7566 0.1748				0.9152 0.2192 0.3721
E localization(mm)				10.000 0 0	22.3607 22.3607 0	17.3205 0 10
	NBI	0.5784 0.1045 0.1378	0.6311 0.0732 0.1189				0.4794 0.1178 0.2574

From top positioning result as can be seen, the E of LPISS _EnergyBe minimum, E _LocalizationOverall performance be better than LORETA and FOCUSS, NBI performance and FOCUSS are suitable, are better than LORETA.These results have shown the sparse source location ability that LIPSS is good, and in that (LORETA in comparison FOCUSS), shows best combination property with current representational method.

2.IOR the source location result of eeg data:

Method: the classical IOR to the Posner design tests the back of making amendment (mainly being the modification that stimulates), adopts amended experiment to stimulate the tested eeg data relevant with IOR that obtain.Eeg data is with EGI system record under the sample rate of 250Hz in 128 roads.Adopt vertex (Cz) electrode as a reference.Behind the record, the processed offline discarded packets contains the data of the excessive electric artefact of eye.15 tested in, the requirements of 13 tested data fit records.Various combination (RR:right cue-location and righttarget-location according to clue and target position of appearing; RL:right cue-location and left target-location; LL:Leftcue-location and left target-location; LR:Left cue-location and righttarget-location), tested record data section (Epoch) is divided into 4 classes: RR, RL, LL and LR.Each data segment length 1.2 seconds originates in the preceding 200ms that target occurs, and lasts till the 1000ms after target occurs.13 effective tested data are carried out population mean according to the combination of clue and target position, obtain 4 classes and bring out current potential: RR, RL, LL and LR.In this work, displaying be that RR is brought out the source location result that current potential carries out.

In IOR experiment, when back 200ms appears in target, a tangible negative wave is arranged by the ERP of target evoked, the composition that is primarily aimed at this moment carries out source location to be tested.In when location, also be the standard real head model that adopts the front, the electrode coordinate of EGI is registrated on the standard header model, adopt BEM to calculate transfer matrix A, its dimension is 128 * 2730, adopts LPISS that 200ms data are constantly carried out source location then.The positioning result of LPISS as shown in Figure 4.

IOR is the brain domain active procedure of a complexity, and the current existing IOR of studies show that is by a lot of brain functioies coefficient result in zone.In utilizing the positioning result of LPISS, (the right frontal eye field in right frontal lobe eye district, RFEF), thalamus (thalamus), right cerebellum (right cerebellum) zone and left occipital lobe (left occipital) position have all navigated to the basic activity zone that stronger activity and current I OR utilize in studying additive method such as fMRI to navigate to and roughly have been consistent.

Claims (4)

1, based on the brain power supply localization method of sparse property constraint, may further comprise the steps:

1) determines transfer matrix A;

2) obtain the EEG signals Y of physical record by multiple tracks EEG signals recording system, carry out pretreatment, determine to carry out the moment that the source is analyzed;

3) initialization source vector x _K-1, k=1, iteration stops the iterations of an error ε and a maximum;

4) upgrade weight diagonal matrix: W _k=(diag (x _K-1));

5) add the constraint of sparse property, utilize optimization method direct estimation q _k: arg min ‖ b-AW _kq _k‖ 2+ λ ‖ q _k‖ _p

6) upgrade source information vector: x _k=W _kq _k

7) stopping criterion for iteration is judged: the source distribution before and after relatively upgrading changes, as ‖ x _k-x _K-1When ‖≤ε or iterations exceeded the maximum iteration time restriction of setting, iteration stopped, x _kBe the final location estimation result in source; Otherwise k=k+1 changes step 2), continue iteration;

Wherein, λ ≠ 0, W _k≠ I.

2, the brain power supply localization method based on the constraint of sparse property as claimed in claim 1 is characterized in that, in the described step 1, comprises step by step following:

A1), the head of object to be measured is carried out MRI or CT scan, obtain the image information of cranial anatomy structure;

A2), extraction step a1) brain part in the image information of gained, then brain is cut apart, extract the functional areas, source (mainly comprising positions such as grey matter, Hippocampus, cerebellum) of brain part again;

A3), with the grid of certain precision with step a2) functional areas, brain source of gained carry out subdivision, determine solution space grid (the locus sequence number that comprises dimension He each grid of solution space);

A4), determine the spatial positional information of each electrode of multiple tracks EEG signals recording system;

A5), determine the model of brain power supply;

A6), utilize step a3) to step a5) in solution space grid, electrode position information and the brain power source model determined, utilize forward modeling method to calculate transfer matrix A, concrete grammar is as follows: the source of placing unit on each solution space position, utilize numerical computation method to calculate the Potential distribution that this unit source produces at the electrode position place, this Potential distribution constitutes the string in the transfer matrix, by that analogy, after all solution space traversals are placed the unit source, just can obtain transfer matrix A.

3, the brain power supply localization method based on the constraint of sparse property as claimed in claim 1 is characterized in that p=1.

4, the brain power supply localization method based on the constraint of sparse property as claimed in claim 1 is characterized in that described step 2) in pretreatment comprise: filtering, baseline correction, eye electricity are removed.

Images