CN105654485A - Magnetic resonance arterial spin labeling sequence partial volume correction method using spatio-temporal information - Google Patents

Abstract

The invention discloses a magnetic resonance arterial spin labeling sequence partial volume correction method using spatio-temporal information and belongs to the medical image partial volume correction technical field. According to the method, based on a situation that the same tested arterial spin labeling sequence and structure image are preprocessed, the four-dimensional data of the arterial spin labeling sequence are constructed; a linear regression (LR) method is utilized to perform three-dimensional spatial correction on the image, and with a correction result adopted as the initial value of the expectation maximization (EM) algorithm, partial volume correction is carried out in a time direction. According to the method (EM-LR) of the invention, the spatial and temporal information of the arterial spin labeling sequence is fully utilized, the problem of the sensitivity of the EM algorithm to an initial value can be solved, the convergence speed of the EM iterative algorithm is increased. As proved by experiments, the invention has excellent edge retention characteristic, and has specificity for small-range perfusion abnormal areas, and at the same time, the method is also applicable to the partial volume correction of other image data with time information.

Description

A kind of magnetic resonance arterial spin labeling Sequence capacity correction method utilizing space time information

Technical field

The invention belongs to the partial volume effect bearing calibration technical field of medical image, particularly to a kind of partial volume correction method of magnetic resonance arterial spin labeling sequence utilizing space time information.

Background technology

Magnetic resonance (MagnaticResonanceImaging, MRI) Perfusion Imaging is used to the blood capillary distribution of reflection tissue and the MRI Examined effect of blood perfusion situation. Arterial spin labeling technology (ArterialSpinLabeling, ASL), by hydrone in magnetization blood as endogenous contrast agent, evaluate cerebral blood perfusion situation, can Non-invasive detection completely, repeatable high, and the cerebral blood flow (CerebralBloodFlow of absolute quantitation can be obtained, CBF), completely not by the impact of blood brain barrier.

ASL technology, using the arterial blood of magnetic marker as endogenous contrast agent, in imaging plane upstream, utilizes the proton of back pulse labelling arterial blood; Postponing after a period of time, blood to be marked enters tissue, and blood and tissue carry out imaging (i.e. labelling image, label) after mass exchange. Label image includes static tissue and the signal of labelling arterial blood. In order to eliminate the signal of static tissue, carrying out additionally once unmarked blood imaging (namely controlling image, control), control image only includes static tissue signal. By label image and control Photographic Subtraction, the error image of gained is only relevant with the labelling arterial blood flowing into imaging plane. The shortcoming of ASL image is in that it is to noise-sensitive. It is, in general, that the difference of label image and control image is about control as the 1��2% of gray value. Therefore, the difference of a pair label/control image can not effectively reflect perfusion situation, it is necessary to about 60 pairs of data obtain rational signal to noise ratio (SignaltoNoiseRatio, SNR).

Meanwhile, as one group of pairing image, label and control image needs to obtain rapidly, continuously. So, fast imaging method (such as EPI method) is widely used in ASL imaging, but the use of this kind of method reduces the spatial resolution of image, makes image be subject to partial volume (PartialVolume, PV) effects serious.At present, (AsllaniI, etal.MagnResonMed, 2008 such as the PV correcting algorithm for ASL image is little, Asllani; 60 (6): 1362-1371.) first by linear regression (LinearRegression, the LR) method of regional area to correct image, but the method makes regional area smooth, and easily causes the CBF error calculated. Afterwards, (LiangX, etal.MagnResonMed, 2013 such as Liang; 69 (2): 531-537.) adopt least trimmed squares method that smoothing effect is carried out post processing. Meanwhile, (ChappellMA, etal.MagnResonMed, 2011 such as Chappell; 65 (4): 1173-1183.) PV of the ASL sequence of multiple reversing times is corrected.

Although ASL data can be carried out a degree of PV correction in spatial domain by said method, but temporal information in ASL sequence being underutilized.

Summary of the invention

In PV bearing calibration for current ASL sequence, the problem that temporal information utilization rate is low, the invention provides a kind of ASL sequence PV bearing calibration utilizing space time information, the method makes full use of the time in ASL sequence and spatial information, sequence is carried out effective PV correction, thus ASL picture quality is greatly improved.

The present invention is achieved through the following technical solutions:

A kind of partial volume correction method of magnetic resonance arterial spin labeling sequence utilizing space time information, comprises the following steps:

(1) MRI data of subjects is gathered, including structure picture and arterial spin labeling sequence;

(2) structure picture and arterial spin labeling sequence are carried out registration;

(3) utilize SPM software that structure picture is split, obtain the probability distribution image of grey matter, white matter and cerebrospinal fluid;

(4) utilize linear regression method that arterial spin labeling sequence is carried out space segment capacity correction, and utilize correction result as the initial value of EM algorithm, time orientation carries out partial volume correction, make full use of the room and time information of arterial spin labeling sequence, it is thus achieved that correct result accurately.

In step (1), the structure picture of collection is T1 or T2 sequence, and arterial spin labeling sequence.

Step (2) is described carries out registration by structure picture and arterial spin labeling sequence, be the MNI coordinate system utilizing SPM to provide as intermediate value to carry out registration.

Step (4) utilizes the 4 D data of arterial spin labeling, along time shaft, the voxel with identical three dimensional space coordinate i is constituted a time arrow; Assuming that all elements is separate in this vector, and grey matter therein and white matter meet Gauss distribution, utilize EM algorithm to estimate line and staff control's central gray and white matter perfusion signal, and concrete method for solving is as follows:

For a certain spatial voxel i, magnetic moment is expressed as:

��M_i=P_iGM��M_iGM+P_iWM��M_iWM(1)

Wherein, P_iGMAnd P_iWMIt it is the probability of mixing voxel central gray and white matter; �� M_iGMWith �� M_iWMRepresenting the perfusion signal of grey matter and white matter, perfusion signal is to be described by the difference DELTA M of label/control image, and the ratios delta M/M of the cerebral blood flow of a certain spatial voxel and �� M and M0 image₀Of equal value; Then the two chamber models according to arterial spin labeling imaging obtain: f_tissue=(�� M/M₀)��F_tissue; Wherein, F_tissueIt is the blood perfusion coefficient relevant to blood brain;

M0 image scans respectively with perfusion image sequence, adopt T1 sequence form, and by with perfusion image registration, convert to and perfusion image formed objects; Therefore assume that M0 image is not affected by the impact of partial volume effect, obtain following relational expression:

$CBF=f_{GM}^P+f_{WM}^P=P_{GM}\left(\frac{{\mathrm{\ΔM}}_{iGM}}{M_0}\right)F_{GM}+P_{WM}\left(\frac{{\mathrm{\ΔM}}_{iWM}}{M_0}\right)F_{WM}---\left(2\right)$

Wherein, P_GMAnd P_WMObtain after carrying out registration by arterial spin labeling sequence and same tested structural images, F_GMAnd F_WMRelevant with the imaging parameters of image;�� M_iGMWith �� M_iWMPerfusion signal for grey matter and white matter; Along time shaft, the voxel with identical three dimensional space coordinate i is constituted a time arrow { Y_it, t=1 ..., T}, T represents the dimension of time arrow; At observed value Y_itIn, the representation in components of grey matter and white matter is:

Y_it=X_itGM+X_itWM(3)

Wherein, X_itGMAnd X_itWMBe respectively with mean forWithVariance isWithStochastic variable, it is assumed that all of T voxel is independent from, then have:

$p\left(X\right|{\stackrel{\‾}{X}}_{iGM},{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iGM}^2,{\mathrm{\σ}}_{iWM}^2)=\underset{t=1}{\overset{T}{\Π}}\left\{p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)\right\}---\left(4\right)$

Perfusion model and formula (4) are combined, obtain:

${\stackrel{\‾}{X}}_{iGM}=P_{iGM}{\mathrm{\ΔM}}_{iGM};{\stackrel{\‾}{X}}_{iWM}=P_{iWM}{\mathrm{\ΔM}}_{iWM};{\mathrm{\σ}}_{iGM}^2=P_{iGM}{S}_{iGM};{\mathrm{\σ}}_{iWM}^2=P_{iWM}{S}_{iWM};$

Wherein, S_iGMAnd S_iWMRepresent grey matter and the variance of white matter perfusion signal respectively;

Assume X_itGMAnd X_itWMGaussian distributed, then formula (4) is transformed to:

$\begin{array}{c}p\left(X\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})=\\ \underset{t=1}{\overset{T}{\Π}}\left\{\left(\frac{1}{\sqrt{2{\mathrm{\πP}}_{iGM}{S}_{iGM}}}{e}^{-\frac{{(X_{itGM}-P_{iGM}{\mathrm{\ΔM}}_{iGM})}^2}{2P_{iGM}{S}_{iGM}}}\right)\left(\frac{1}{\sqrt{2{\mathrm{\πP}}_{iWM}{S}_{iWM}}}{e}^{-\frac{{(X_{itWM}-P_{iWM}{\mathrm{\ΔM}}_{iWM})}^2}{2P_{iWM}{S}_{iWM}}}\right)\right\}\end{array}---\left(5\right)$

In EM algorithm, the observed value Y of t voxel_itIt it is an incomplete stochastic variable; X_itGMAnd X_itWMWhat represent is line and staff control's information complete in the t voxel, is a complete variable;

It is being integrated for condition with formula (1), is setting up imperfect variable { Y_itAnd complete variable { X_itGMAnd { X_itWMBetween the relation of probability distribution, such as following formula:

$\begin{array}{c}p\left(Y_{it}\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})\\ ={\∫}_{\{Y_{it}=X_{itGM}+X_{itWM}\}}p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)dX\end{array}---\left(6\right)$

Wherein, i represents a certain voxel three-dimensional space position, { Y_itRepresent and there is the voxel of same spatial location i constituted a time arrow; { X_itGMAnd { X_itWMRespectively observed value { Y_itThe component of central gray and white matter.

The optimal solution obtaining step of formula (6) is as follows:

Adopting EM algorithm to ask for the greatest hope of complete Variable Conditions probability distribution, E-step is that M-step is used for asking for expectation maximum for the log-likelihood of variable is estimated:

E-step: design conditions probability expectation p (X | ��), wherein

Conditional expectation is expressed as:

$\begin{array}{c}Q\left(\mathrm{\Θ}\right|{\mathrm{\Θ}}^{\left(n\right)})=E_{Y_{it}=X_{itGM}+X_{itWM}}\[\ln \left(p\right(X|\mathrm{\Θ}\left)\right)|Y,{\mathrm{\Θ}}^{\left(n\right)}\]\\ E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[X_{itWM}^2-2P_{iWM}{\mathrm{\ΔM}}_{iGM}{X}_{itGM}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]+\\ E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[X_{itWM}^2-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{X}_{itWM}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]\\ =-\frac{1}{2}\underset{t}{\Σ}\left\{\begin{array}{c}\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]+\\ \ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\end{array}\right\}\end{array}---\left(7\right)$

Derivation according to conditional expectation, draws:

$\begin{array}{c}{X}_{itGM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGMt}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+\frac{P_{iGM}{S}_{iGM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(8\right)$

$\begin{array}{c}{X}_{itWM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+\frac{P_{iWM}{S}_{iWM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(9\right)$

${\left(X_{itGM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itGM}^2|Y_{it},\mathrm{\Θ}\left(n\right)\]={\left(X_{itGM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(10\right)$

${\left(X_{itWM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itWM}^2|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]={\left(X_{itWM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(11\right)$

M-step: make conditional probability expected value maximize by n+1 iteration, the average in line and staff control's model is by maximizing n+1 iteration of conditional probabilityObtain:

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iGM}}{|}_{{\mathrm{\ΔM}}_{iGM}={\mathrm{\ΔM}}_{iGM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itGM}^{\left(n\right)}}{T\·P_{iGM}}---\left(12\right)$

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iWM}}{|}_{{\mathrm{\ΔM}}_{iWM}={\mathrm{\ΔM}}_{iWM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iWM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itWM}^{\left(n\right)}}{T\·P_{iWM}}---\left(13\right)$

$S_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itGM}^2\right)}^{\left(n\right)}-2X_{itGM}^{\left(n\right)}{P}_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}\right)}^2\]}{T\·P_{iGM}}---\left(14\right)$

$S_{iWM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itWM}^2\right)}^{\left(n\right)}-2X_{itWM}^{\left(n\right)}{P}_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}\right)}^2\]}{T\·P_{iWM}}---\left(15\right).$

Step (4) concrete operations are as follows:

1) for every width dimensional perfusion image, linear regression method is used to carry out space segment capacity correction;

2) based on the space segment capacity correction result of linear regression method, the voxel with same spatial location i is constituted time arrow;

3) for time arrow, the line and staff control that greatest hope method carries out on time shaft is utilized to estimate;

4) to locus i+1, the 2nd is repeated) and the 3rd) step, until view picture perfusion image is all corrected.

Compared with prior art, the present invention has following useful technique effect:

The ASL Sequence capacity correction method based on space time information of the present invention, first, gathers the MRI data of subjects, the structure picture of collection and ASL sequence is carried out registration; Then, utilize SPM software that structure picture is split, generate grey matter (graymatter respectively, GM), white matter (whitematter, and cerebrospinal fluid (cerebrospinalfluid WM), CSF) probability distribution image, it is possible to the probability distribution information of GM and WM is provided for follow-up PV correction, is conducive to estimating the parameter of line and staff control; Again, utilize LR method that ASL sequence carries out space PV correction; Finally, utilizing the free-air correction result initial value as greatest hope (expectationmaximization, EM) algorithm of LR method, and be corrected on time orientation, thus accurately estimating cerebral blood flow value, contributing to follow-up data analysis. The method (EM-LR) takes full advantage of the room and time information of arterial spin labeling sequence, compensate for EM algorithm to initial value sensitive issue, accelerates the convergence rate of EM iterative algorithm.The partial volume correction method utilizing space time information provided by the invention, is applicable not only to ASL sequence, is also applied for other partial volume correction with temporal information image data.

Accompanying drawing explanation

Fig. 1 is the flow chart of steps that the present invention utilizes the magnetic resonance arterial spin labeling Sequence capacity correction method of space time information;

Fig. 2 be use do not correct, analog data is corrected the intermediate layer result figure of the GM cerebral blood flow scattergram obtained by EM method, LR method and EM-LR method;

Fig. 3 is under different noise conditions, the result that the CBF area-of-interest (ROI) of GM is analyzed; Wherein, (a) noise intensity is that ROI when 2.5 analyzes, and (b) noise intensity is that ROI when 6.5 analyzes;

Fig. 4 is EM, LR and EM-LR method, respectively in the correction result at HT and Low perfusion regional center line; Wherein, a () and (b) is respectively when noise intensity is 2.5, the CBF change of GM on the centrage in HT and Low perfusion region, c () and (d), respectively when noise intensity is 6.5, on the centrage in HT and Low perfusion region, the CBF of GM changes;

Fig. 5 is for the correction result of analog simulation 2, LR, EM and EM-LR.

Detailed description of the invention

Below in conjunction with specific embodiment, the present invention is described in further detail, and the explanation of the invention is not limited.

Referring to Fig. 1, a kind of magnetic resonance arterial spin labeling Sequence capacity correction method utilizing space time information disclosed by the invention, comprise the steps:

(1) after carrying out tested preparation, tested MRI data is gathered, including structure picture and ASL sequence;

(2) structure picture and ASL sequence are carried out registration;

(3) utilize SPM (StatisticalParametricMapping) software that structure picture is split, obtain the probability distribution image of GM, WM and CSF;

(4) utilize LR method that ASL sequence carries out space PV correction, and it can be used as the initial value of EM algorithm, time orientation carries out PV correction (EM-LR algorithm), thus making full use of the room and time information of ASL sequence, it is thus achieved that correct result accurately.

In said method, the structure picture of step (1) described collection can be T1 or T2 sequence. Method according to arterial blood inversion marks is different, ASL technology can be broadly divided into two types: continuous way arterial spin labeling (ContinuousArterySpinLabeling, and pulsed arterial spin labeling (PulsedArterySpinLabeling, PASL) CASL). The present invention can carry out PV correction for any type of ASL sequence, and generally, the multiple scanning data of ASL sequence are approximately 60 width, in order to obtain rational signal to noise ratio.

Structure picture and ASL sequence described in step (2) carry out registration, utilize SPM MNI (MontrealNeurologicalInstitute) coordinate system provided as intermediate value to carry out registration. When using SPM to read in structure picture and ASL data, producing the transition matrix between image coordinate system and MNI coordinate system, what MNI coordinate was corresponding is standard form, it is possible to for the registration of structure picture and ASL data.

Structure picture segmentation described in step (3), utilizes the normal process that SPM software provides to split, generates the probability distribution image of GM, WM and CSF respectively. The probability distribution image that segmentation obtains, can provide the probability distribution information of GM and WM, be conducive to estimating the parameter of line and staff control for follow-up PV correction.

The time orientation PV based on EM algorithm described in step (4) corrects.Utilize the 4 D data of ASL, along time shaft, the voxel with identical three dimensional space coordinate i be may be constructed a time arrow. Assuming that all elements is separate in this vector, and GM and WM therein meets Gauss distribution, and in line and staff control, GM and WM irrigates signal to utilize EM algorithm to estimate, detailed process is as follows:

Owing to the space rate respectively of ASL image is low, mix, at one, the combined effect may including GM, WM and CSF in voxel. It is generally believed that the CBF in ASL sequence is not contributed by CSF, therefore, for a certain spatial voxel i, magnetic moment can be expressed as

��M_i=P_iGM��M_iGM+P_iWM��M_iWM(1)

Wherein, P_iGMAnd P_iWMIt it is the probability of GM and WM in mixing voxel; �� M_iGMWith �� M_iWMRepresent the signal of GM and WM.

By ASL image image-forming principle it can be seen that perfusion signal is to be described by the difference DELTA M of label/control image, and ratio (the �� M/M of the CBF of a certain voxel and �� M and M0 image₀) of equal value. Therefore, can obtain according to two chamber models of ASL imaging, f_tissue=(�� M/M₀)��F_tissue, wherein F_tissueIt is the blood perfusion coefficient relevant to blood brain. M0 image scans respectively with perfusion image sequence, what adopt is the form of T1 sequence, and by with perfusion image registration, convert to and perfusion image formed objects, therefore M0 image is only small by the impact of PV effect, therefore this research hypothesis M0 image is not affected by the impact of PV effect. Thus, it is possible to obtain following relational expression:

$CBF=f_{GM}^P+f_{WM}^P=P_{GM}\left(\frac{{\mathrm{\ΔM}}_{iGM}}{M_0}\right)F_{GM}+P_{WM}\left(\frac{{\mathrm{\ΔM}}_{iWM}}{M_0}\right)F_{WM}---\left(2\right)$

In above-mentioned cerebral blood flow measure equation, P_GMAnd P_WM(step 3) can be obtained, F after carrying out registration by ASL sequence and same tested structural images_GMAnd F_WMRelevant with the imaging parameters of image, next it is how to calculate the perfusion signal delta M of GM and WM_iGMWith �� M_iWM��

In order to obtain �� M_iGMWith �� M_iWM, utilize the 4 D data of ASL, along time shaft, the voxel with identical three dimensional space coordinate i constituted a time arrow { Y_it, t=1 ..., T}, T represents the dimension of time arrow. At observed value Y_itIn, the component of GM and WM is represented by X_itGMAnd X_itWM, therefore,

Y_it=X_itGM+X_itWM(3)

Wherein, X_itGMAnd X_itWMBe respectively with mean forWithVariance isWithStochastic variable. Assume what all of T voxel was independent from, then

$p\left(X\right|{\stackrel{\‾}{X}}_{iGM},{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iGM}^2,{\mathrm{\σ}}_{iWM}^2)=\underset{t=1}{\overset{T}{\Π}}\left\{p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)\right\}---\left(4\right)$

Perfusion model and above-mentioned formula are combined, it is possible to obtain Wherein S_iGMAnd S_iWMRepresent the variance of GM and WM signal respectively.

Assume X_itGMAnd X_itWMGaussian distributed, then formula (4) can transform to:

$\begin{array}{c}p\left(X\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})=\\ \underset{t=1}{\overset{T}{\Π}}\left\{\left(\frac{1}{\sqrt{2{\mathrm{\πP}}_{iGM}{S}_{iGM}}}{e}^{-\frac{{(X_{itGM}-P_{iGM}{\mathrm{\ΔM}}_{iGM})}^2}{2P_{iGM}{S}_{iGM}}}\right)\left(\frac{1}{\sqrt{2{\mathrm{\πP}}_{iWM}{S}_{iWM}}}{e}^{-\frac{{(X_{itWM}-P_{iWM}{\mathrm{\ΔM}}_{iWM})}^2}{2P_{iWM}{S}_{iWM}}}\right)\right\}\end{array}---\left(5\right)$

In EM algorithm, the observed value Y of t voxel_itIt it is an incomplete stochastic variable. And X_itGMAnd X_itWMWhat represent is line and staff control's information complete in the t voxel, is therefore a complete variable. It is being integrated with formula (3) for condition, it is possible to represent imperfect variable { Y_itAnd complete variable { X_itGMAnd { X_itWMBetween the relation of probability distribution,

$\begin{array}{c}p\left(Y_{it}\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})\\ ={\∫}_{\{Y_{it}=X_{itGM}+X_{itWM}\}}p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)dX\end{array}---\left(6\right)$

In order to obtain the optimal solution of equation (6), EM algorithm is used to ask for the greatest hope of complete Variable Conditions probability distribution, and wherein E-step is that M-step is used for asking for expectation maximum for the log-likelihood of variable is estimated.

E-step: design conditions probability expectation p (X | ��), whereinConditional expectation can be expressed as:

$\begin{array}{c}Q\left(\mathrm{\Θ}\right|{\mathrm{\Θ}}^{\left(n\right)})=E_{Y_{it}=X_{itGM}+X_{itWM}}\[\ln \left(p\right(X|\mathrm{\Θ}\left)\right)|Y,{\mathrm{\Θ}}^{\left(n\right)}\]\\ E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[X_{itGM}^2-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{X}_{itGM}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]+\\ E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[X_{itWM}^2-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{X}_{itWM}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]\\ =-\frac{1}{2}\underset{t}{\Σ}\left\{\begin{array}{c}\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]+\\ \ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\end{array}\right\}\end{array}---\left(7\right)$

Derivation according to conditional expectation, it can be deduced that:

$\begin{array}{c}{X}_{itGM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGMt}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+\frac{P_{iGM}{S}_{iGM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(8\right)$

$\begin{array}{c}{X}_{itWM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+\frac{P_{iWM}{S}_{iWM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(9\right)$

${\left(X_{itGM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itGMt}^2|Y_{it},\mathrm{\Θ}\left(n\right)\]={\left(X_{itGM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(10\right)$

${\left(X_{itWM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itWM}^2|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]={\left(X_{itWM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(11\right)$

M-step: make conditional probability expected value maximize by n+1 iteration. Average in line and staff control's model can pass through to maximize n+1 iteration of conditional probabilityObtain, namely

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iGM}}{|}_{{\mathrm{\ΔM}}_{iGM}={\mathrm{\ΔM}}_{iGM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itGM}^{\left(n\right)}}{T\·P_{iGM}}---\left(12\right)$

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iWM}}{|}_{{\mathrm{\ΔM}}_{iWM}={\mathrm{\ΔM}}_{iWM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iWM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itWM}^{\left(n\right)}}{T\·P_{iWM}}---\left(13\right)$

$S_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itGM}^2\right)}^{\left(n\right)}-2X_{itGM}^{\left(n\right)}{P}_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}\right)}^2\]}{T\·P_{iGM}}---\left(14\right)$

$S_{iWM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itWM}^2\right)}^{\left(n\right)}-2X_{itWM}^{\left(n\right)}{P}_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}\right)}^2\]}{T\·P_{iWM}}---\left(15\right)$

Pass through said method, it is possible to use the temporal information of ASL sequence, obtain the perfusion signal delta M of GM and the WM of position i_iGMWith �� M_iWM, but owing to EM algorithm is more sensitive to initial value, use relatively accurate initial value, can the convergence of accelerating algorithm, improve the calculating accuracy of EM algorithm.

For this, the present invention proposes two kinds of initial value plans of establishment:

1) EM method: utilize the spatial prior information found in document as the initial value of algorithm, namely during E-step, adopt ${\mathrm{\ΔM}}_{iGM}^{\left(0\right)}=50,S_{iGM}^{\left(0\right)}=1,{\mathrm{\ΔM}}_{iWM}^{\left(0\right)}=10,S_{iWM}^{\left(0\right)}=1,$ EM method is initialized, and utilizes the line and staff control that EM method (formula 8-15) carries out on time shaft to estimate.

2) EM-LR method: will combine based on the free-air correction of LR method and the time adjustment based on EM algorithm, utilizes the correction result of LR method to provide initial value for EM algorithm. Detailed process is as follows:

1. for every width dimensional perfusion image, LR method (AsllaniI, etal.MagnResonMed, 2008 are used; 60 (6): 1362-1371.) space PV correction is carried out;

2. based on the free-air correction result of LR method, the voxel with same spatial location i constitutes time arrow;

3. for time arrow, the line and staff control that EM method (formula 8-15) carries out on time shaft is utilized to estimate;

4. to locus i+1, the 2nd is repeated) and the 3rd) step, until view picture perfusion image is all corrected.

Utilizing the ASL sequence PV of space time information to correct as it has been described above, just complete, the method takes full advantage of time and the spatial information of ASL sequence, compensate for EM algorithm to initial value sensitive issue, accelerates the convergence rate of EM iterative algorithm.

Analog simulation 1:

Utilizing analog simulation to generate height and irrigate region, using method and traditional method disclosed in this patent that model is carried out PV correction, thus evaluating the effectiveness of this method.

Phantom generation step:

On the basis using SPM software that magnetic resonance T1 sequence is standardized and to split, obtain the probability distribution image of GM, WM and CSF. ASL sequence is simulated in the following way:

1) GM and WM image being standardized, image is sized to 60 �� 72 �� 60, and voxel size is 3 �� 3 �� 3mm³;

2) the CBF value of WM is set to �� M=20ml/100g/min;

3) the CBF value of GM is set to �� M=60ml/100g/min, and chooses the bulbous region that 2 radiuses are 5 voxels, is respectively set to Low perfusion and HT district, and is decided to be 30 and 90ml/100g/min;

4) by formula (1), view picture dimensional perfusion image is obtained;

5) two kinds of Gaussian noise intensity 2.5 (value conventional in document) and 6.5 (noises that in document, the ASL sequence of report is the strongest) are adopted to be added in dimensional perfusion imaging;

6) the multiple scanning number of times (60 times) according to conventional ASL sequence, repeats 1-5 step, generates two groups of four-dimension ASL sequences.

Simulation result:

With reference to Fig. 2, use do not correct, analog data is corrected by EM method, LR method and EM-LR method. It can be seen that EM, LR and EM-LR method is better than uncorrected result, especially at the intersection of GM and CSF, but the result that LR method obtains has obvious smoothing effect at the edge in high/low perfusion region.

With reference to Fig. 3, area-of-interest (RegionofInterest, ROI) is analyzed, and wherein, (a) noise intensity is that ROI when 2.5 analyzes, and (b) noise intensity is that ROI when 6.5 analyzes. For the CBF measurement accuracy of Quantitative Comparison distinct methods, GM average CBF concordance in GM probability distribution is adopted to be evaluated. Based on the probability distribution image of GM, choose 9 regions, i.e. P_GMAt [10%��20%], [20%��30%] ... [90%��100%] scope, and calculate the average CBF value of GM in each scope respectively. Fig. 3 illustrates under different noise level, and the ROI of analog data analyzes result.Similar with the observed result of Fig. 2, result that EM, LR and EM-LR method obtains is close with actual value, but in the relatively low region of GM probability distribution, the result that EM-LR method and LR method obtain is better.

With reference to Fig. 4, in order to study the effect in edge maintenance of EM, LR, EM-LR method, choose two sections through center at HT and Low perfusion region respectively. Wherein, a () and (b) is respectively when noise intensity is 2.5, the CBF change of GM on the centrage in HT and Low perfusion region, c () and (d), respectively when noise intensity is 6.5, on the centrage in HT and Low perfusion region, the CBF of GM changes. Fig. 4 illustrates and intersects on hatching line at section with the intermediate layer shown in Fig. 2, the comparison between the estimated value of actual value and EM, LR and EM-LR method. It can be seen that EM and EM-LR method has good edge retention performance, it is possible to accurately estimate cerebral blood flow value while effectively keeping details.

Analog simulation 2:

In order to study the inventive method, zonule and/or slight perfusion are changed the calibration accuracy in region, by " analog simulation 1 " the 3rd) step is revised as: the CBF value of GM is set to �� M=60ml/100g/min, and wherein, 1. choose the bulbous region that radius is 5 voxels, �� M=65ml/100g/min (slight HT) is set; 2. select the cubical area of 3 �� 3 �� 3, �� M=55ml/100g/min (zonule, slight Low perfusion) is set; 3. select the cubical area of 2 �� 2 �� 2, �� M=65ml/100g/min (zonule, and slight HT) is set.

Simulation result:

With reference to Fig. 5, although the region studied is smaller, and perfusion changes little, but uses EM method and EM-LR method can obtain and well correct result, and LR method is difficult to little perfusion abnormality region is corrected.

By above-mentioned two model emulation, EM disclosed by the invention (use experience value is as the initial value of EM algorithm) and EM-LR (uses the free-air correction result initial value as EM algorithm of LR method) can carry out PV correction to ASL data, and has specificity for the perfusion abnormality that region is less. But due to the individual variation of people, empirical value may differ greatly with actual value, so utilizing LR method more reasonable as the initial value (EM-LR) of EM algorithm.

By simulation results show, the partial volume correction method utilizing space time information provided by the invention, it is possible to arterial spin labeling sequence is carried out effective partial volume correction, the perfusion abnormality region less simultaneously for scope has specificity. Meanwhile, the method is also applied for the partial volume correction of other image datas with temporal information.

Claims (6)

1. the partial volume correction method of the magnetic resonance arterial spin labeling sequence utilizing space time information, it is characterised in that comprise the following steps:

(1) MR data of subjects is gathered, including structure picture and arterial spin labeling sequence;

(2) structure picture and arterial spin labeling sequence are carried out registration;

(3) utilize SPM software that structure picture is split, obtain the probability distribution image of grey matter, white matter and cerebrospinal fluid;

(4) utilize linear regression method that arterial spin labeling sequence is carried out space segment capacity correction, and utilize correction result as the initial value of EM algorithm, time orientation carries out partial volume correction, make full use of the room and time information of arterial spin labeling sequence, it is thus achieved that correct result accurately.

2. the partial volume correction method of the magnetic resonance arterial spin labeling sequence utilizing space time information according to claim 1, it is characterised in that in step (1), the structure picture of collection is T1 or T2 sequence, and arterial spin labeling sequence.

3. the partial volume correction method of the magnetic resonance arterial spin labeling sequence utilizing space time information according to claim 1, it is characterized in that, step (2) is described carries out registration by structure picture and arterial spin labeling sequence, be the MNI coordinate system utilizing SPM to provide as intermediate value to carry out registration.

4. the partial volume correction method of the magnetic resonance arterial spin labeling sequence utilizing space time information according to claim 3, it is characterized in that, step (4) utilizes the 4 D data of arterial spin labeling, along time shaft, the voxel with identical three dimensional space coordinate i is constituted a time arrow; Assuming that all elements is separate in this vector, and grey matter therein and white matter meet Gauss distribution, utilize EM algorithm to estimate line and staff control's central gray and white matter perfusion signal, and concrete method for solving is as follows:

For a certain spatial voxel i, magnetic moment is expressed as:

��M_i=P_iGM��M_iGM+P_iWM��M_iWM(1)

Wherein, P_iGMAnd P_iWMIt it is the probability of mixing voxel central gray and white matter; �� M_iGMWith �� M_iWMRepresenting the perfusion signal of grey matter and white matter, perfusion signal is to be described by the difference DELTA M of label/control image, and the ratios delta M/M of the cerebral blood flow of a certain spatial voxel and �� M and M0 image₀Of equal value; The two chamber models according to arterial spin labeling imaging obtain: f_tissue=(�� M/M₀)��F_tissue;

Wherein, F_tissueIt is the blood perfusion coefficient relevant to blood brain;

M0 image scans respectively with perfusion image sequence, adopt T1 sequence form, and by with perfusion image registration, convert to and perfusion image formed objects; Therefore assume that M0 image is not affected by the impact of partial volume effect, obtain following relational expression:

$CBF=f_{GM}^P+f_{WM}^P=P_{GM}\left(\frac{{\mathrm{\ΔM}}_{iGM}}{M_0}\right)F_{GM}+P_{WM}\left(\frac{{\mathrm{\ΔM}}_{iWM}}{M_0}\right)F_{WM}---\left(2\right)$

Wherein, P_GMAnd P_WMObtain after carrying out registration by arterial spin labeling sequence and same tested structural images, F_GMAnd F_WMRelevant with the imaging parameters of image; �� M_iGMWith �� M_iWMPerfusion signal for grey matter and white matter; Along time shaft, the voxel with identical three dimensional space coordinate i is constituted a time arrow { Y_it, t=1 ..., T}, T represents the dimension of time arrow; At observed value Y_itIn, the representation in components of grey matter and white matter is:

Y_it=X_itGM+X_itWM(3)

Wherein, X_itGMAnd X_itWMBe respectively with mean forWithVariance isWithStochastic variable, it is assumed that all of T voxel is independent from, then have:

$p\left(X\right|{\stackrel{\‾}{X}}_{iGM},{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iGM}^2,{\mathrm{\σ}}_{iWM}^2)=\underset{t=1}{\overset{T}{\Π}}\left\{p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)\right\}---\left(4\right)$

Perfusion model and formula (4) are combined, obtain:

${\stackrel{\‾}{X}}_{iGM}=P_{iGM}{\mathrm{\ΔM}}_{iGM};{\stackrel{\‾}{X}}_{iWM}=P_{iWM}{\mathrm{\ΔM}}_{iWM};{\mathrm{\σ}}_{iGM}^2=P_{iGM}{S}_{iGM};{\mathrm{\σ}}_{iWM}^2=P_{iWM}{S}_{iWM};$

Wherein, S_iGMAnd S_iWMRepresent grey matter and the variance of white matter perfusion signal respectively;

Assume X_itGMAnd X_itWMGaussian distributed, then formula (4) is transformed to:

$\begin{array}{c}p\left(X\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})=\\ \underset{t=1}{\overset{T}{\Π}}\left\{\left(\frac{1}{\sqrt{2\mathrm{\π}}{P}_{iGM}{S}_{iGM}}{e}^{-\frac{{(X_{itGM}-P_{iGM}{\mathrm{\ΔM}}_{iGM})}^2}{2P_{iGM}{S}_{iGM}}}\right)\left(\frac{1}{\sqrt{2\mathrm{\π}}{P}_{iWM}{S}_{iWM}}{e}^{-\frac{{(X_{itWM}-P_{iWM}{\mathrm{\ΔM}}_{iWM})}^2}{2P_{iWM}{S}_{iWM}}}\right)\right\}\end{array}---\left(5\right)$

In EM algorithm, the observed value Y of t voxel_itIt it is an incomplete stochastic variable; X_itGMAnd X_itWMWhat represent is line and staff control's information complete in the t voxel, is a complete variable;

It is being integrated for condition with formula (1), is setting up imperfect variable { Y_itAnd complete variable { X_itGMAnd { X_itWMBetween the relation of probability distribution, such as following formula:

$\begin{array}{c}p\left(Y_{it}\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})=\\ ={\∫}_{(Y_{it}=X_{itGM}+X_{itWM})}p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)dX\end{array}---\left(6\right)$

Wherein, i represents a certain voxel three-dimensional space position, { Y_itRepresent and there is the voxel of same spatial location i constituted a time arrow; { X_itGMAnd { X_itWMRespectively observed value { Y_itThe component of central gray and white matter.

5. the partial volume correction method of the magnetic resonance arterial spin labeling sequence utilizing space time information according to claim 4, it is characterised in that the optimal solution obtaining step of formula (6) is as follows:

Adopting EM algorithm to ask for the greatest hope of complete Variable Conditions probability distribution, E-step is that M-step is used for asking for expectation maximum for the log-likelihood of variable is estimated:

E-step: design conditions probability expectation p (X | ��), wherein

Conditional expectation is expressed as:

$\begin{array}{c}Q\left(\mathrm{\Θ}\right|{\mathrm{\Θ}}^{\left(n\right)})=E_{Y_{it}=X_{itGM}+X_{itWM}}\[\ln \left(p\right(X|\mathrm{\Θ}\left)\right)|Y,{\mathrm{\Θ}}^{\left(n\right)}\]\\ =E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[X_{itGM}^2-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{X}_{itGM}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]+\\ E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[X_{itWM}^2-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{X}_{itWM}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]\\ =-\frac{1}{2}\underset{t}{\Σ}\left\{\begin{array}{c}\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]+\\ \ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\end{array}\right\}\end{array}---\left(7\right)$

Derivation according to conditional expectation, draws:

$\begin{array}{c}{X}_{itGM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGMt}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+\frac{P_{iGM}{S}_{iGM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(8\right)$

$\begin{array}{c}{X}_{itWM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+\frac{P_{iWM}{S}_{iWM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(9\right)$

${\left(X_{itGM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itGMt}^2|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]={\left(X_{itGM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(10\right)$

${\left(X_{itWM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itGM}^2|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]={\left(X_{itGM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(11\right)$

M-step: make conditional probability expected value maximize by n+1 iteration, the average in line and staff control's model is by maximizing n+1 iteration of conditional probabilityObtain:

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iGM}}{|}_{{\mathrm{\ΔM}}_{iGM}={\mathrm{\ΔM}}_{iGM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itGM}^{\left(n\right)}}{T\·P_{iGM}}---\left(12\right)$

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iWM}}{|}_{{\mathrm{\ΔM}}_{iWM}={\mathrm{\ΔM}}_{iWM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iWM}^{(n+1)}=\frac{{\Σ}_{t=1}^T{X}_{itWM}^{\left(n\right)}}{T\·P_{iWM}}---\left(13\right)$

$S_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itGM}^2\right)}^{\left(n\right)}-2X_{itGM}^{\left(n\right)}{P}_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}\right)}^2\]}{T\·P_{iGM}}---\left(14\right)$

$S_{iWM}^{(n+1)}=\frac{{\Σ}_{t=1}^T\[{\left(X_{itWM}^2\right)}^{\left(n\right)}-2X_{itWM}^{\left(n\right)}{P}_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}\right)}^2\]}{T\·P_{iWN}}---\left(15\right).$

6. the partial volume correction method of the magnetic resonance arterial spin labeling sequence utilizing space time information according to claim 4, it is characterised in that step (4) concrete operations are as follows:

1) for every width dimensional perfusion image, linear regression method is used to carry out space segment capacity correction;

2) based on the space segment capacity correction result of linear regression method, the voxel with same spatial location i is constituted time arrow;

3) for time arrow, the line and staff control that greatest hope method carries out on time shaft is utilized to estimate;

4) to locus i+1, the 2nd is repeated) and the 3rd) step, until view picture perfusion image is all corrected.