CN108470335B - Multi-correlation source scanning imaging method based on brain source space segmentation - Google Patents

Multi-correlation source scanning imaging method based on brain source space segmentation Download PDF

Info

Publication number
CN108470335B
CN108470335B CN201810239440.0A CN201810239440A CN108470335B CN 108470335 B CN108470335 B CN 108470335B CN 201810239440 A CN201810239440 A CN 201810239440A CN 108470335 B CN108470335 B CN 108470335B
Authority
CN
China
Prior art keywords
source
activity
brain
virtual
distribution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810239440.0A
Other languages
Chinese (zh)
Other versions
CN108470335A (en
Inventor
陈丹
蔡畅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN201810239440.0A priority Critical patent/CN108470335B/en
Publication of CN108470335A publication Critical patent/CN108470335A/en
Application granted granted Critical
Publication of CN108470335B publication Critical patent/CN108470335B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30016—Brain

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

The invention discloses a multi-correlation source scanning imaging method based on brain source space segmentation, which aims at the imaging analysis of data acquired by a human brain at a sleeping moment, adopts a brand-new distributed brain source reconstruction scanning method-a multi-layer region combined source reconstruction SOCHU method, and overcomes the defect that the analysis cannot be carried out by an empirical Bayes estimation and convex plane boundary distributed tomography method (namely, a Champagne method). The problem that the correlation source cannot be reconstructed by adaptive beamforming is solved, and the constraint that methods such as Champagne need to depend on data before stimulation is avoided. The SOCHU method adopts the idea of scanning one by one with the self-adaptive beam forming method to solve the source activity, when any potential active source is solved, other incoming signals are not required to be suppressed, and the other incoming signals are estimated by using uninteresting virtual sources, so that the solution of the potential source activity is ensured, and the influence of correlation is removed.

Description

Multi-correlation source scanning imaging method based on brain source space segmentation
Technical Field
The invention belongs to the interdisciplinary field of brain science and information technology, relates to an imaging method for brain source space segmentation, and particularly relates to a distributed multi-correlation source scanning imaging method based on brain source space segmentation.
Background
Magnetoencephalography (MEG) and electroencephalography (EEG) are two popular methods for non-destructive examination of brain activity by measuring magnetic and electric fields on the scalp surface generated by recording brain activity, and then performing imaging studies on brain activity. The process of human brain activity may be achieved by brain-derived activity reconstruction of the acquired MEG and EEG data. And combining the brain source activity model, the brain structure information and the information of the sensor system, and solving a leadfield guide matrix through a Maxwell equation set so as to describe the relationship between the real brain source activity and the data acquired by the sensor. The invention provides a brand-new distributed brain Source reconstruction scanning method, namely a multilayer regional joint Source reconstruction (SOCHU) method, aiming at the problem that data acquired by a human brain at a sleeping moment can not be analyzed by an empirical Bayes estimation and convex plane boundary distributed tomography method (namely a Champagne method).
Disclosure of Invention
Aiming at the problems, the invention provides a brain Source space segmentation-based distributed multi-correlation Source scanning imaging method, namely a multi-layer region joint Source reconstruction (SOCHU) method, aiming at solving the problem that a correlation Source cannot be reconstructed by self-adaptive beam forming and avoiding the constraint that the method such as Champagne and the like depends on data before stimulation.
The technical scheme adopted by the invention is as follows: a multi-correlation source scanning imaging method based on brain source space segmentation comprises the following steps:
step 1: generating a multilayer region joint source reconstruction probability model;
step 2: estimating the distribution of the multi-layer region combined source hyper-parameters;
and step 3: determining a multi-layer region joint source virtual source;
and 4, step 4: the multi-layer region is associated with a source-residual minimum variance adaptive beamforming method.
The concept of scanning one by one to solve the source activity by the SOCHU and the adaptive beam forming method is different from the latter in that: when any potential active source is solved, other incoming signals are not required to be suppressed, and the other incoming signals are estimated by using uninteresting virtual sources, so that the solution of the potential source activity is ensured, and the influence of correlation is removed. The SOCHU iterative cost function is also sampled based on the Bayes theory and the idea of convex function boundary, and the convergence speed (relative to the EM method) is ensured. Based on different assumptions about virtual source activity, the present invention proposes three methods of SOCHU: (1) assuming that all virtual source activities have the same distribution, the same as the Minimum Norm (MN), called SOCHU _ MN; (2) the virtual source is divided into different areas, and the different areas have different source activity distributions, which are called SOCHU _ tiles; (3) the virtual sources are identified by singular value decomposition and each virtual source has a different source activity distribution, called SOCHU _ SVD.
Compared with the prior art, the invention has the beneficial effects that: aiming at the problem that data collected in sleep and the like cannot be analyzed by a Champagne method, the invention provides a brand-new distributed brain Source reconstruction scanning method, namely a multi-layer regional joint Source reconstruction (SOCHU) method, which can be used for solving the problem that a correlation Source cannot be reconstructed by self-adaptive beam forming, avoids the problem that the Champagne method and the like need to depend on the constraint of data before stimulation, and has certain effect.
Drawings
FIG. 1 is a flow chart of an embodiment of the present invention;
FIG. 2 is a graph illustrating the convergence comparison of two types of iterative methods when estimating the source activity at location i according to an embodiment of the present invention.
Detailed Description
In order to facilitate the understanding and implementation of the present invention for those of ordinary skill in the art, the present invention is further described in detail with reference to the accompanying drawings and examples, it is to be understood that the embodiments described herein are merely illustrative and explanatory of the present invention and are not restrictive thereof.
Referring to fig. 1, the multi-correlation source scanning imaging method based on brain source space segmentation provided by the present invention includes the following steps:
step 1, generating a universal probability model by the SOCHU method. The method comprises the following specific steps:
step 1.1: assume that the sensor acquires data for an extracerebral MEG (or EEG) at time t, denoted as y (t) ═ y1(t),y2(t),…,yM(t)]T(i ═ 1,2, …, M), where yi(t) represents the extracerebral data acquired by the sensor at time t, and M is the number of channels of the acquisition point. By discrete segmentation of the brain activity region, while assuming each voxel is potentially active as a brain source. The brain-derived activity at time t is
Figure BDA0001604797580000021
N is the number of potential brain source activities of the segmentation, and the source activity s of the ith element at the time ti(t) is represented by dcVector of dimensions
Figure BDA0001604797580000022
Where d iscIndicating the number of source activity directions.
Then when solving for the brain source activity scan at location i, assuming that the entire brain activity is summed up by the activity of the ith element and other virtual sources, the model is:
y(t)＝lisi(t)+g\iv\i(t)+ε (1)
here, the
Figure BDA0001604797580000031
A leadfield matrix for position i, liThe k-th column of (a) represents the source activity at position i in k-direction, the extra-brain sensor acquires the intensity of its activity-producing signal. In general, when source reconstructing MEG data, assume dcLet d be assumed when source reconstructing EEG data 2cIt solves in 2D and 3D space corresponding to the source activity direction. The scholars put forward a plurality of solutions to the leadfield matrix l based on the Maxwell equation setiThe method of (1). g\i＝[g1,g2,,gQ-1]Is the leadfield matrix corresponding to the rest Q-1 virtual source activities, which comprise all brain activities except the i source, background interference and noise, Q is the sum of the number of the virtual source and the real source i activities, v\i(t)＝[v1(t),v2(t),…,vQ-1(t)]TIs the activity intensity at time t corresponding to the virtual source, the jth virtual source activity
Figure BDA0001604797580000032
Finally, the parameters ε are independent of each other and the distribution obeys a normal distribution
Figure BDA0001604797580000033
(ii) a gaussian noise signal;
step 1.2: all virtual sources are decomposed into R regions according to brain anatomy or functional structure, assuming region j contains pjThe virtual sources, the partitioning used here is a non-coincident partitioning, i.e. each virtual source belongs to only one region. Equation (1) can be expressed as:
Figure BDA0001604797580000034
gj,mand vj,m(t) respectively correspond to the virtual sources numbered m located in the region j,
Figure BDA0001604797580000035
representing the summation of the extracerebral data generated by all virtual sources belonging to region j. If all position variables are assigned probability distributions, the model is a probability model.
Step 1.3: the noise epsilon of the formula (2) is assumed to be known, because the position noise covariance matrix can be used for estimating noise information (estimated by data before stimulation or task data, namely, noise) by performing Variational Bayesian Factor Analysis (VBFA) on the acquired signals to solve residual signals of the acquired signals. It is then assumed that both the source activity parameter and the noise activity are gaussian distributed. At the same time, it is assumed that they are independent of each other and have the same random probability distribution at any time.
Let it be assumed here that the a priori covariance distribution of the source activity i is dc×dcIs matrix phiiHere, the subscript i is used to indicate that this scan is the source activity s for location iiThe solution is carried out, so that:
Figure BDA0001604797580000041
virtual source v\i(t) the prior distribution of activity is assumed to be:
Figure BDA0001604797580000042
here psijIs dc×dcRepresenting the covariance matrix ψ of virtual sources having the same distribution located in region jj. For simple calculation, if the ith virtual source is located in region j, we define the ith virtual source covariance matrix as γi＝ψjTherefore, the following can be obtained:
Figure BDA0001604797580000043
here:
Figure BDA0001604797580000044
if the brain-derived activity at location i is considered as a region with only one brain-derived activity i, a simpler and more uniform general model can be obtained by simplifying equation (2) as follows:
Figure BDA0001604797580000045
wherein:
Figure BDA0001604797580000046
here, it should be emphasized that in the model of equation (7), only the active sources in the first region are the true source activities solved by the current scan, so f1＝liAnd x1＝siThe key of the calculation is, and the rest source activities are pseudo sources.
Step 1.4: suppose y (t)k) May be represented as y (k), x (t)k) Can be expressed as X (k), and the brain source activity for the entire time window is expressed as X ═ X (t) >1),x(t2),…,x(tK)]The data collected by the sensor outside the brain is represented as Y ═ Y (t)1),y(t2),…,y(tK)]Assume that its prior distribution is:
Figure BDA0001604797580000051
where the prior covariance matrix Ω is dcQ×dcBlock diagonal matrix for Q:
Figure BDA0001604797580000052
here, ΩiIs dc×dcWhen i is 1, Ωi＝φ1The source activity prior covariance matrix corresponding to position i, Ω when i is 2,3, …, Qi＝γi-1The a priori covariance matrix corresponding to the i-1 th virtual source activity. Since it has been assumed that the distribution of the noise ε is known, as
Figure BDA0001604797580000053
Therefore, the conditional probability distribution p (Y | X) can be expressed as:
Figure BDA0001604797580000054
to estimate the distribution of the source X, it is necessary to solve the posterior probability distribution p (X | Y) of the model, defined as:
Figure BDA0001604797580000055
the precision and mean of the posterior probability obtained by Bayesian estimation in combination with equation (9), equation (11) and equation (12) are:
Figure BDA0001604797580000056
Figure BDA0001604797580000057
the posterior mean can also be expressed as:
Figure BDA0001604797580000058
wherein:
∑y＝∑ε+fΩfT (16)
for the covariance matrix of the model data, the result of equation (15) can be expressed as:
Figure BDA0001604797580000059
the updated formula for j source activity of equation (15) is:
Figure BDA00016047975800000510
equation (18) is the time series update rule for the SOCHU method.
And 2, estimating the hyperparameter distribution in the SOCHU method. The method comprises the following specific steps:
step 2.1: brain-derived activity parameter xkThe estimate of (c) can be updated by either equation (14) or equation (18), but there still exists an unknown hyperparameter Ω. In bayesian estimation, the estimation of the hyperparameter Ω is performed by maximizing the marginal probability distribution p (Y | Ω), which is expressed as:
Figure BDA0001604797580000061
although the value of the hyperparameter Ω can be optimally solved by maximizing the marginal probability distribution p (Y | Ω), the second term log | Σ on the right of equation (19)yThe presence of | makes it difficult to maximize equation (19), a process commonly referred to as class 2 maximum likelihood estimation or empirical bayesian estimation. In fact, the solution to this problem can also be estimated by the expectation maximization (EM method), which uses the M step of the EM method:
Figure BDA0001604797580000062
step 2.2: since the convergence speed of the EM method is slow, the cost function of equation (19) is modified here to accelerate by using the same auxiliary function (convex function boundary) as the Champagne method. Due to log | ∑yI is a convex function of a hyperparameter omega, so d is introducedc×dcAuxiliary parameter matrix Λj(j ═ 1,2, …, Q), satisfying:
Figure BDA0001604797580000063
here Λ0For scalar values, then defining a secondary cost function
Figure BDA0001604797580000064
Comprises the following steps:
Figure BDA0001604797580000065
it always satisfies:
Figure BDA0001604797580000071
so if the value of the over-parameter omega enables the auxiliary cost function
Figure BDA0001604797580000072
The value of (2) is increased, and the hyperparameter omega also increases the value of the marginal probability logp (Y | omega), so the problem is converted into an auxiliary cost function
Figure BDA0001604797580000073
The optimization of the maximum value comprises the following steps:
Figure BDA0001604797580000074
the updating rule of the hyper-parameter is derived by an auxiliary cost function:
Figure BDA0001604797580000075
setting equation (25) to 0, one obtains:
Figure BDA0001604797580000076
since the requirement of Ω satisfies equation (23), the result is:
Figure BDA0001604797580000077
equation (27) is the update rule of the hyper-parameter Ω.
Step 2.3: finding the auxiliary parameter ΛjDue to formulas
Figure BDA0001604797580000078
For convex function log | ∑yA strict ceiling function of |, sojAlways being hyperplane function log | ∑yThe tangent plane slope of l. Therefore, ΛjCan be updated by applying log | ∑ toySolving the slope of the hyperplane to obtain:
Figure BDA0001604797580000079
here omegajFor block diagonal matrix elements that exceed the parameter Ω, the following needs to solve for the average covariance matrix for each partition region, since the covariance matrix for region j is given by the assumption that all elements within each region have the same covariance matrix:
Figure BDA00016047975800000710
then, the values of the formula (18), the formula (27), the formula (28) and the formula (29) solved this time are used as initial values of the next iteration loop, and the iteration operation is continued, so that the convergence of the auxiliary function (or the edge distribution) is known.
Step 2.4: the value of the edge distribution function is increased at each iteration, and the value monitoring is carried out by formula (19) or a strict lower limit auxiliary function formula (22) of the edge distribution likelihood function. Because the likelihood function depends on the brain source activity of the position i, when any position activity source is subjected to iterative calculation, the corresponding likelihood function is used for monitoring the iterative process of the method, and in addition, each time the calculation of one element is finished, the time sequence of the position activity source can be solved through a formula (18). The complexity of each iteration operation is a linear function of Q, so the method has moderate calculation complexity. The convergence speed of the auxiliary cost function based on the convex function boundary is much faster than that of the cost function based on the EM method, fig. 2(b) is an iteration convergence comparison graph of the two methods when the elements at the position i are solved, and it can be seen from the graph that the auxiliary function based on the convex function boundary completes convergence in 25 iterations, and the EM method still does not completely converge after more than 150 iterations.
In summary, the hyperparameter Ω is obtained by iterative operations of formula (18), formula (27), formula (28) and formula (29)jIs determined, which ensures that the value of the edge distribution likelihood function logp (Y | Ω) is increased or kept constant during each iteration, when it converges, i.e. the source activity covariance matrix Ω of the location ij. And then, carrying out the same iterative solution on elements at other positions in the source space until the estimation of all element covariance matrixes is completed, and finally solving the element activity energy and the time sequence of all the positions in the source space through a formula (17). Therefore, this brain source reconstruction method becomes SOCHU.
And step 3: determination of virtual sources in the SOCHU method. The method comprises the following specific steps:
step 3.1: based on the formula (6) model and Bayesian estimation, the distribution of all the source activities and the corresponding activity time sequence are obtained by iteratively solving all the potential activity sources segmented in the source activity space one by one, however, the settings of the virtual source by the SOCHU method are unknown all the time, and the effective virtual source and the corresponding leadfield matrix g thereof will be performed below\iAnd (4) solving. When evaluating the source activity at location i, it is assumed that the virtual source and its leadfield matrix are the same as the source activity at non-location i and its leadfield matrix, i.e., g\i＝l\iHere l\i＝[l1,…,li-1,li+1,…,lN]N is the number of elements into which the source space is divided, with Q ═ N. Under the assumption, if all virtual sources are set to have the same distribution and the same covariance matrix, that is, the virtual sources are not subjected to region segmentation R ═ 1, all virtual sources belong to the same region, and a model based on a multi-layer region joint source based on the model is called SOCHU _ MN; suppose that the virtual source is divided into a number of regions R>The virtual sources belonging to different areas have different distributions and covariance matrixes, the virtual sources in the same area have the same distribution, and the multi-layer area combined source set based on the distribution is called SOCHU _ tiles.
However, in many cases there is not enough a priori knowledge to simply or efficiently partition the virtual source into different regions, another approach is to perform a singular value decomposition on the leadfield matrix corresponding to elements other than position i, as follows:
l\i＝BΔD* (30)
assuming that the number of virtual sources is Q-1, the leadfield matrix g corresponding to the virtual sources\iIs g\iBased on this assumption, different virtual sources correspond to different gaussian distributions, B (: 1: Q-1) Δ (1: Q-1 ). In fact, the optimal number Q-1 of virtual sources can be determined according to the geometric condition of the element, such as a line segment segmentation method. The SOCHU based on this Singular Value Decomposition (SVD) is referred to as SOCHU _ SVD.
And 4, step 4: the relationship of the SOCHU and the minimum variance adaptive beam forming method is solved. The method comprises the following specific steps:
step 4.1: minimum Variance Adaptive Beamforming (MVAB) is currently one of the most widely used methods in brain source activity reconstruction. MVAB method by
Figure BDA0001604797580000091
The activity of a dipole source at a location i at a time t, where WiThe spatial domain filter for element i is defined as:
Figure BDA0001604797580000092
as is known, MVAB cannot locate more than two strongly correlated source activities, and although many scholars attempt to solve this problem, it is only applicable for certain special cases. The MVAB method can also be regarded as a special form of the SOCHU method, and the MVAB is the SOCHU method calculated in one iteration. Model y (t) l of formula (1)isi(t)+g\iv\i(t) + ε, it can be concluded that the sum of interference and noise is z when the source activity estimation is performed on the element at position in＝g\iv\i(t) + ε, the covariance matrix of the sampled signals y (t) may be expressed as ∑y＝∑ε+fΩfT. Suppose that the model parameters have been estimated as ∑εAnd phiiThen the SOCHU has a source activity time sequence of position i elements of
Figure BDA0001604797580000093
Wherein:
Figure BDA0001604797580000094
here phiiFor the prior distribution value of the location i source activity, equation (32) and equation (31) are the same after transformation.
Step 4.2: assuming that the covariance matrix of the sampled data satisfies the condition of infinite data sampling points
Figure BDA0001604797580000095
Using matrix inversion theorem, RyyCan be expressed as
Figure BDA0001604797580000096
Figure BDA0001604797580000097
Thus, there are:
Figure BDA0001604797580000098
the estimate of Θ is used below as
Figure BDA0001604797580000099
It assumes that the source activity prior distribution variance at location i is approximately equal to 0. Therefore, the method comprises the following steps:
Figure BDA0001604797580000101
the last step uses formula (33).
In specific implementation, those skilled in the art can implement automatic operation of the above processes by using computer software technology.
It should be understood that parts of the specification not set forth in detail are well within the prior art.
It should be understood that the above description of the preferred embodiments is given for clarity and not for any purpose of limitation, and that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (3)

1. A multi-correlation source scanning imaging method based on brain source space segmentation is characterized by comprising the following steps:
step 1: generating a multilayer region joint source reconstruction probability model;
the specific implementation of the step 1 comprises the following substeps:
step 1.1: assume that the sensor acquires data for an extracerebral MEG or EEG at time t, denoted as y (t) ═ y1(t)，y2(t)，…，yM(t)]TWherein y isi(t) represents extracerebral data acquired by the sensor at time t, i is 1, 2. By discrete segmentation of brain active regions, while assuming each voxel as a potential brain sourceThe activity of the brain source at time t is
Figure FDA0003322729740000011
N is the number of potential brain source activities of the segmentation, and the source activity s of the ith element at the time ti(t) is represented by dcVector of dimensions
Figure FDA0003322729740000012
Where d iscRepresenting the number of the source moving directions;
then when solving for the brain source activity scan at location i, assuming that the entire brain activity is summed up by the activity of the ith element and other virtual sources, the model is:
y(t)＝lisi(t)+g\iv\i(t)+ε (1)
here, the
Figure FDA0003322729740000013
A leadfield matrix for position i, liWhen the source at the position i is in unit source activity in the k direction, the sensor outside the brain acquires the intensity of the signal generated by the activity; g\i＝[g1，g2，...，gQ-1]Is the leadfield matrix corresponding to the rest Q-1 virtual source activities, which comprise all brain activities except the i source, background interference and noise, Q is the sum of the number of the virtual source and the real source i activities, v\i(t)＝[v1(t)，v2(t)，…，vQ-1(t)]TIs the activity intensity at time t corresponding to the virtual source, the jth virtual source activity
Figure FDA0003322729740000014
The parameters epsilon are independent of each other and the distribution follows a normal distribution
Figure FDA0003322729740000015
(ii) a gaussian noise signal;
step 1.2: all virtual sources are assigned according to the brain anatomy or functional structureDividing the image into R regions without overlap, and assuming that the region j includes pjA virtual source; equation (1) is then expressed as:
Figure FDA0003322729740000016
gj，mand vj，m(t) respectively correspond to the virtual sources numbered m located in the region j,
Figure FDA0003322729740000017
representing the summation of the extracerebral data generated by all virtual sources belonging to the region j; if all the position variables are endowed with probability distribution, the model is a probability model;
step 1.3: assuming that the noise epsilon of the formula (2) is known, assuming that the source activity parameter and the noise activity are Gaussian distributions, they are independent of each other and have the same random probability distribution at any time;
let the prior covariance distribution of the source activity i be dc×dcIs matrix phiiThe subscript i indicates that this scan is a source activity s for location iiThe solution is carried out, so that:
Figure FDA0003322729740000021
virtual source v\i(t) the prior distribution of activity is assumed to be:
Figure FDA0003322729740000022
here psijIs dc×dcRepresenting the covariance matrix ψ of virtual sources having the same distribution located in region jj；
If the ith virtual source is located in the region j, defining the covariance matrix of the ith virtual source as gammai＝ψjTherefore, the following components are obtained:
Figure FDA0003322729740000023
here:
Figure FDA0003322729740000024
if the brain-derived activity at location i is considered as a region with only one brain-derived activity i, the general model for equation (2) is simplified as follows:
Figure FDA0003322729740000025
wherein:
Figure FDA0003322729740000026
in the model of equation (7), only the active sources in the first region are the real source activity for the current scan solution, so f1＝liAnd x1＝siThe key of the calculation is, and the other source activities are pseudo sources;
step 1.4: suppose y (t)k) Expressed as y (k), x (t)k) Denoted X (k), the brain source activity for the entire time window is denoted X ═ X (t)1)，x(t2)，…，x(tK)]The data collected by the sensor outside the brain is represented as Y ═ Y (t)1)，y(t2)，…，y(tK)]Assume that its prior distribution is:
Figure FDA0003322729740000031
where the prior covariance matrix Ω is dcQ×dcBlock diagonal matrix for Q:
Figure FDA0003322729740000032
here, ΩiIs dc×dcWhen i is 1, Ωi＝φ1The source activity prior covariance matrix corresponding to position i, Ω when i is 2,3, …, Qi＝γi-1A prior covariance matrix corresponding to the i-1 th virtual source activity; since it has been assumed that the distribution of the noise ε is known, as
Figure FDA0003322729740000033
Therefore, the conditional probability distribution p (Y | X) is expressed as:
Figure FDA0003322729740000034
to estimate the distribution of the source X, it is necessary to solve the posterior probability distribution p (X | Y) of the model, defined as:
Figure FDA0003322729740000035
the precision and the mean value of the posterior probability obtained by Bayesian estimation and combined with the formula (9), the formula (11) and the formula (12) are as follows:
Figure FDA0003322729740000036
Figure FDA0003322729740000037
the posterior means are expressed as:
Figure FDA0003322729740000038
wherein:
∑y＝∑ε+fΩfT (16)
for the covariance matrix of the model data, the result of equation (15) is expressed as:
Figure FDA0003322729740000041
the updated formula for j source activity of equation (15) is:
Figure FDA0003322729740000042
formula (18) is a time series updating rule of the multi-layer region joint source reconstruction probability model;
step 2: estimating the distribution of the multi-layer region combined source hyper-parameters;
the specific implementation of the step 2 comprises the following substeps:
step 2.1: in bayesian estimation, the estimation of the hyperparameter Ω is performed by maximizing the marginal probability distribution p (Y | Ω), which is:
Figure FDA0003322729740000043
step 2.2: modifying the cost function of the formula (19), and accelerating by using an auxiliary function which is the same as the Champagne method; due to log | ∑yI is a convex function of a hyperparameter omega, so d is introducedc×dcAuxiliary parameter matrix ΛjJ is 1,2, …, Q, satisfying:
Figure FDA0003322729740000044
here Λ0For scalar values, then defining a secondary cost function
Figure FDA0003322729740000045
Comprises the following steps:
Figure FDA0003322729740000046
it always satisfies:
Figure FDA0003322729740000047
so if the value of the over-parameter omega enables the auxiliary cost function
Figure FDA0003322729740000048
The value of (2) is increased, and the hyperparameter omega also increases the value of the marginal probability logp (Y | omega), so the problem is converted into an auxiliary cost function
Figure FDA0003322729740000049
The optimization of the maximum value comprises the following steps:
Figure FDA0003322729740000051
the updating rule of the hyper-parameter is derived by an auxiliary cost function:
Figure FDA0003322729740000052
setting equation (25) to 0, yields:
Figure FDA0003322729740000053
since the requirement of Ω satisfies equation (23), the result is:
Figure FDA0003322729740000054
the formula (27) is an updating rule of the hyper-parameter omega;
step 2.3: finding the auxiliary parameter ΛjThe update rule of (1);
Λjby applying the update rule of (1) to log | ∑ySolving the slope of the hyperplane to obtain:
Figure FDA0003322729740000055
here omegajFor block diagonal matrix elements that exceed the parameter Ω, the following needs to solve for the average covariance matrix for each partition region, since the covariance matrix for region j is given by the assumption that all elements within each region have the same covariance matrix:
Figure FDA0003322729740000056
then taking the values of the formula (18), the formula (27), the formula (28) and the formula (29) solved at this time as initial values of the next iteration loop, continuing the iteration operation, and knowing the convergence of the auxiliary function or edge distribution;
and step 3: determining a multi-layer region joint source virtual source;
the specific implementation process of the step 3 is as follows:
when evaluating the source activity at location i, it is assumed that the virtual source and its leadfield matrix are the same as the source activity at non-location i and its leadfield matrix, i.e., g\i＝l\iHere l\i＝[l1，…，li-1，li+1，…，lN]N is the number of elements into which the source space is divided, with Q ═ N; under this assumption, if all virtual sources are set to have the same distribution and the same covariance matrix, i.e. they are not subjected to region segmentation R ═ 1, all virtual sources belong to the same regionThe domain, the model based on the multi-layer region joint source is called SOCHU _ MN; if the number R of the virtual source partition areas is larger than 1, the virtual sources belonging to different areas have different distributions and covariance matrixes, the virtual sources in the same area have the same distribution, and the multi-layer area combined source set based on the distribution is called SOCHU _ tiles;
and 4, step 4: a multi-layer region combined source residue minimum variance adaptive beam forming method is connected;
the specific implementation of the step 4 comprises the following substeps:
step 4.1: minimum variance adaptive beamforming MVAB method
Figure FDA0003322729740000061
The activity of a dipole source at a location i at a time t, where WiThe spatial domain filter for element i is defined as:
Figure FDA0003322729740000062
model y (t) l of formula (1)isi(t)+g\iv\i(t) + ε, the sum of interference and noise is z when the source activity estimation is performed for the element at position in＝g\iv\i(t) + ε, the covariance matrix of the sampled signals y (t) may be expressed as ∑y＝∑ε+fΩfT(ii) a Suppose that the model parameters have been estimated as ∑εAnd phiiThen the source activity time series of the multi-layer regional federated source to the location i element is
Figure FDA0003322729740000063
Wherein:
Figure FDA0003322729740000064
here phiiThe prior distribution value of the source activity of the position i is changed by the formula (32) and the formula (31)The same is obtained after the replacement;
step 4.2: assuming that the covariance matrix of the sampled data satisfies the condition of infinite data sampling points
Figure FDA0003322729740000065
Using matrix inversion theorem, RyyCan be expressed as
Figure FDA0003322729740000066
Figure FDA0003322729740000067
Therefore, there are:
Figure FDA0003322729740000068
the estimate of Θ is used below as
Figure FDA0003322729740000069
It is assumed that the source activity prior distribution variance at location i is approximately equal to 0; therefore, the method comprises the following steps:
Figure FDA00033227297400000610
Figure FDA0003322729740000071
the last step uses formula (33).
2. The multi-correlation source scanning imaging method based on brain source space segmentation according to claim 1, characterized in that in step 2.1, the estimation of the value of the hyper-parameter Ω is performed by the expectation maximization EM method, with the M step of the EM method:
Figure FDA0003322729740000072
3. the multi-correlation source scanning imaging method based on brain source space segmentation according to claim 1, characterized in that: if there is not enough a priori knowledge to simply or efficiently partition the virtual source into different regions, the leadfield matrix corresponding to elements other than position i is subjected to a singular value decomposition as follows:
l\i＝BΔD* (30)
assuming that the number of virtual sources is Q-1, the leadfield matrix g corresponding to the virtual sources\iIs g\iB (: 1: Q-1) Δ (1: Q-1), based on which the different virtual sources correspond to different gaussian distributions; the optimal number of virtual sources Q-1 is determined according to the geometric conditions of the elements, and the multi-layer region joint source decomposed based on the singular values is called SOCHU _ SVD.
CN201810239440.0A 2018-03-22 2018-03-22 Multi-correlation source scanning imaging method based on brain source space segmentation Active CN108470335B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810239440.0A CN108470335B (en) 2018-03-22 2018-03-22 Multi-correlation source scanning imaging method based on brain source space segmentation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810239440.0A CN108470335B (en) 2018-03-22 2018-03-22 Multi-correlation source scanning imaging method based on brain source space segmentation

Publications (2)

Publication Number Publication Date
CN108470335A CN108470335A (en) 2018-08-31
CN108470335B true CN108470335B (en) 2022-02-15

Family

ID=63264618

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810239440.0A Active CN108470335B (en) 2018-03-22 2018-03-22 Multi-correlation source scanning imaging method based on brain source space segmentation

Country Status (1)

Country Link
CN (1) CN108470335B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109965869B (en) * 2018-12-16 2021-09-10 北京工业大学 MI-EEG identification method based on brain source domain space
CN111317466A (en) * 2019-07-03 2020-06-23 重庆邮电大学 Electroencephalogram signal imaging method and system and computer equipment
CN113995422B (en) * 2021-10-15 2022-11-25 苏州大学 Transient brain power source positioning method and system based on nonnegative block sparse Bayesian learning

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7801591B1 (en) * 2000-05-30 2010-09-21 Vladimir Shusterman Digital healthcare information management
CN105249964A (en) * 2015-11-10 2016-01-20 东南大学 Multimodal brain function reconstruction assessment method based on magnetoencephalogram and diffusion tensor imaging

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101531994B1 (en) * 2013-09-04 2015-06-29 한국과학기술연구원 Apparatus and method for selectively collecting electroencephalogram data through motion recognition
GB201409766D0 (en) * 2014-06-02 2014-07-16 Cambridge Entpr Ltd Signal processing methods

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7801591B1 (en) * 2000-05-30 2010-09-21 Vladimir Shusterman Digital healthcare information management
CN105249964A (en) * 2015-11-10 2016-01-20 东南大学 Multimodal brain function reconstruction assessment method based on magnetoencephalogram and diffusion tensor imaging

Also Published As

Publication number Publication date
CN108470335A (en) 2018-08-31

Similar Documents

Publication Publication Date Title
Armanious et al. Unsupervised medical image translation using cycle-MedGAN
Yoon et al. Motion adaptive patch-based low-rank approach for compressed sensing cardiac cine MRI
CN108470335B (en) Multi-correlation source scanning imaging method based on brain source space segmentation
Adler et al. Learning to solve inverse problems using Wasserstein loss
US11969283B2 (en) Generating synthetic electron density images from magnetic resonance images
CN112690775B (en) Bayes-based imaging system for focal zone with abnormal brain activity of children
CN109003232B (en) Medical MRI image denoising method based on frequency domain scale smoothing Shearlet
US10102649B2 (en) Iterative image subset processing for image reconstruction
Du et al. Robust statistical recognition and reconstruction scheme based on hierarchical Bayesian learning of HRR radar target signal
Kim et al. Iterative approach of dual regression with a sparse prior enhances the performance of independent component analysis for group functional magnetic resonance imaging (fMRI) data
US10303849B2 (en) One gate reconstruction
CN108416822B (en) Bayesian estimation-based multi-level and multi-scale tomography method
CN109584322B (en) Shearlet medical PET image denoising method based on frequency domain direction smoothing
US11741643B2 (en) Reconstruction of dynamic scenes based on differences between collected view and synthesized view
CN115048963A (en) System and method for simultaneously solving brain source activity and noise
López et al. Single MEG/EEG source reconstruction with multiple sparse priors and variable patches
US10217249B2 (en) Image reconstruction using approximate message passing methods
Le Folgoc et al. Sparse bayesian registration
CN117751388A (en) Method of noninvasive medical tomography with uncertainty estimation
Yan et al. MAP image reconstruction using intensity and line processes for emission tomography data
Subudhi et al. Context dependent fuzzy associated statistical model for intensity inhomogeneity correction from magnetic resonance images
Swetha et al. Sparse feature aware noise removal technique for brain multiple sclerosis lesions using magnetic resonance imaging
López et al. Random location of multiple sparse priors for solving the MEG/EEG inverse problem
Ouzir et al. Cardiac motion estimation in ultrasound images using spatial and sparse regularizations
Mejia et al. A spatial template independent component analysis model for subject-level brain network estimation and inference

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant